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Abstract 

The global nonlinear behavior of finite discretizations for constant time steps and fixed 
or adaptive grid spacings is studied using tools from dynamical systems theory. Detailed 
analysis of commonly used temporal and spatial discretizations for simple model problems 
is presented. The role of dynamics in the understanding of long time behavior of numerical 
integration and the nonlinear stability, convergence, and reliability of using time-marching 
approaches for obtaining steady-state numerical solutions in computational fluid dynamics 
(CFD) is explored. The study is complemented with examples of spurious behavior observed 
in steady and unsteady CFD computations. The CFD examples were chosen to illustrate 
non-apparent spurious behavior that was difficult to detect without extensive grid and temporal 
refinement studies and some knowledge from dynamical systems theory. Studies revealed the 
various possible dangers of misinterpreting numerical simulation of realistic complex flows that 
are constrained by available computing power. In large scale computations where the physics 
of the problem under study is not well understood and numerical simulations are the only viable 
means of solution, extreme care must be taken in both computation and interpretation of the 
numerical data. The goal of this paper is to explore the important role that dynamical systems 
theory can play in the understanding of the global nonlinear behavior of numerical algorithms 
and to aid the identification of the sources of numerical uncertainties in CFD. 


l An earlier version of this paper was published as an internal report N ASA Technical Memorandum 
110398 April 1990. 
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1. Introduction 

This paper is an outgrowth of the NASA Technical Memorandum 110398, April 1996, 
entitled “Nonlinear Dynamics and Numerical Uncertainties in CFD.” This expanded version 
includes approximately 30% new material. Many sections have been rewritten and many 
sections have been shortened or deleted to accommodate more practical examples in spurious 
behavior (numerical artifacts) of unsteady computational fluid dynamics (CFD) simulations. 
The new examples presented in Sections 5.4 - 5.6 were chosen based on their non-apparent 
numerical uncertainties that were difficult to detect without extensive spatial and temporal 
refinement studies and some knowledge from dynamical systems. 

The authors’ views and experience in the application of nonlinear dynamical systems 
theory to improve the understanding of global nonlinear behavior of finite discretizations and 
their connection to numerical uncertainties in CFD are reviewed. Simple nonlinear model 
equations are used to illustrate how the recent advances in nonlinear dynamical systems theory 
can provide new insights and further the understanding of nonlinear effects on the asymptotic 
behavior of numerical algorithms commonly used in CFD. The discussion is complemented with 
CFD examples containing spurious behavior in steady and unsteady flows. Although this paper 
is intended primarily for computational fluid dynamicists, it can be useful for computational 
scientists, physicists, engineers and computer scientists who have a need for reliable numerical 
simulation. 

Since the late 1980’s, many CFD related journals imposed an editorial policy statement 
on numerical uncertainty which pertained mainly to the accuracy issue. However, the study 
of numerical uncertainties in practical computational physics encompasses very broad subject 
areas. These include but are not limited to (a) problem formulation and modeling, (b) type, 
order of accuracy, nonlinear stability, and convergence of finite discretizations, (c) limits 
and barriers of existing finite discretizations for highly nonlinear stiff problems with source 
terms and forcing, and/or for wave propagation phenomena, (d) numerical boundary condition 
procedures, (e) finite representation of infinite domains (f) solution strategies in solving 
the nonlinear discretized equations, (g) procedures for obtaining the steady-state numerical 
solutions, (h) grid quality and grid adaptations, (i) multigrids, and (j) domain decomposition 
(zonal or multicomponent approach) in solving large problems. See, forexample, Mehta (1995), 
Melnik et al. (1994), Cosner et al. (1995), Demuren & Wilson (1994), Marvin (1993), Marvin 
& Holst (1990) and references cited therein on guidelines for code verification, validation and 
certification. At present, some of the numerical uncertainties can be explained and minimized by 
traditional numerical analysis and standard CFD practices. However, such practices might not 
be sufficient for strongly nonlinear and/or stiff problems. Examples of this type of problem are 
combustion, direct numerical simulations, high speed and reacting flows, and certain turbulence 
models in Navier-Stokes computations. We believe that a good understanding of the nonlinear 
behavior of numerical schemes being used should be an integral part of code verification and 
validation. See Jackson (1989) for the definition of highly (or genuinely) nonlinear problems. 

A major stumbling block in genuinely nonlinear studies is that unlike the linear model 
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equations used for conventional stability and accuracy considerations in time-dependent partial 
differential equations (PDEs), there is no equivalent unique nonlinear model equation for 
nonlinear hyperbolic and parabolic PDEs for fluid dynamics. A numerical method behaving in 
a certain way for a particular nonlinear differential equation (DE) (PDE or ordinary differential 
equation (ODE)) might exhibit a different behavior for a different nonlinear DE even though the 
DEs are of the same type. On the other hand, even for simple nonlinear model DEs with known 
solutions, the discretized counterparts can be extremely complex, depending on the numerical 
methods. Except in special cases, there is no general theory at the present time to characterize 
the various nonlinear behaviors of the underlying discretized counterparts. Most often, the only 
recourse is a numerical approach. Under this constraint, whenever analytical analysis of the 
discretized counterparts is not possible, the associated dynamics of numerics such as bifurcation 
phenomena and asymptotic behavior are obtained numerically using supercomputers. It is 
hoped that we can encourage numerical analysts to construct practical algorithms (to avoid 
spurious dynamics) based on the numerical phenomena observed using supercomputers to 
balance advances of computations and analyses. We also hope that it will strengthen the 
interface of numerical analysis with practical CFD applications and motivate CFD researchers 
who are looking for new approaches and solutions to new or old but challenging problems. 

The term “discretized counterparts” is used to mean the finite difference equations 
resulting from finite discretizations of the underlying DEs. Here “dynamics” is used loosely 
to mean the dynamical behavior of nonlinear dynamical systems (continuum or discrete) 
and “numerics” is used loosely to mean the numerical methods and procedures in solving 
dynamical systems. We emphasize here that in the study of the dynamics of numerics, unless 
otherwise stated, we always assume the continuum (governing equations) is nonlinear. 

Outline : A rather detailed background, motivation and subtleties of the subject will be 
given in Section II due to the relatively new yet interdisciplinary nature of this research topic 
for CFD. Particular attention is paid to the isolation of the different nonlinear behavior and 
spurious dynamics due to some of the numerical uncertainties that are observed in practical 
CFD computations. The background material includes the connection between (continuum) 
nonlinear dynamical systems and fluid dynamics, between nonlinear discrete dynamical systems 
and CFD, and between nonlinear dynamics and time-marching approaches. With the aid of 
elementary examples. Section III reviews the fundamentals of spurious behaviors of explicit 
and implicit time discretizations and spatial discretizations. Sections IV and V illustrate typical 
CFD computations that exhibit spurious behavior similar to that in the elementary examples. 
The majority of the material in these two sections have been reported in Yee & Sweby (1996a) 
and Yee et al. (1997). Sections 4.2.1, 5.1, 5.2, 5.5 and 5.6 were written by the original 
contributors of the respective work. Section 5.4 is the joint work of the first author with John 
Torczynski of the Sandia National Laboratories (Yee et al. 1997). The discussion is divided into 
transient and steady-state computations with several examples for each category. Section IV is 
mainly concerned with convergence rate and spurious behavior of time-marching to the steady 
states of high-resolution shock-capturing methods. Section V is concerned with momentum 
spikes and post- shock oscillations in slowly moving shocks, “numerically induced chaotic 
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transient” computations, and spurious behavior due to underresolved grid and semi- implicit 
temporal discretizations. The concluding remarks section discusses how to use existing tools 
in bifurcation theory to minimize convergence to the wrong steady state or asymptote. These 
tools are based on the combined knowledge of recent advances in the dynamics of the physical 
equations and the dynamics of the underlying finite discretizations. 


n. Background & Motivation 

Since the late 1970’s there have been important developments and breakthroughs in the 
theory of nonlinear dynamical systems. There was also an explosion of journal and conference 
papers, texts and reference books on the subject in the 1980’s and early 1990’s. See for example, 
Guckenheimer & Holmes (1983), Thompson & Stewart (1988), Seydel (1988), Hales & Kocak 
(1991), Stewart (1990), Wiggins (1990), and Hoppensteadt (1993). During the early 1980’s, 
a new area of applied mathematics emerged from the interaction of dynamical system theory 
and numerical analysis. These developments addressed mainly mathematical principles and 
their applications of numerics in the understanding of the dynamics of DEs without discussing 
the connection between dynamics and numerics for initial value problems (IVPs) and initial 
boundary value problems (IB VPs). There was, however, some discussion on this connection for 
boundary value problems (BVPs) (Doedel & Beyn 1981, Beyn & Lorenz 1982, Doedel 1986, 
Beyn 1987, Shubin et al. 1981, Kellogg et al. 1980, Peitgen et al. 1981 and Shreiber & Keller 
1983). Studies of BVPs of the elliptic type continue to the present day. See, for example, the last 
four SIAM Conferences on Dynamical Systems (1990, 1992, 1995, 1997) and the Proceedings 
of IMA Conference on Dynamics of Numerics and Numerics of Dynamics, July 31 - Aug. 2, 
1990, Bristol, England. 

Why is it important to understand the connection between dynamics and numerics for BVPs, 
IVPs and IB VPs? It stems from the fact that it is a standard practice to use numerics to discover 
dynamical properties of continuous systems. As a matter of fact, much of what we know about 
specific dynamical systems is usually obtained from numerical experiments. One not only can 
visualize the dynamics and the bifurcation phenomena associated with numerics, but in most 
of the cases the dynamical behavior of the DEs is not amenable otherwise. Consequently, 
developments concerned with the connection between dynamics and numerics are necessary 
to bridge the gap in a better understanding of the dynamics of numerics and the numerics of 
dynamics. 

In the late 1980’s, developments concerned with the connection between dynamics and 
numerics for IVPs and IB VPs slowly emerged. See for example, Mitchell & Griffiths (1985), 
Griffiths & Mitchell (1988), Pruffer (1985), Iserles (1988), Lorenz (1989), Sanz-Sema (1985, 
1990) and Sanz-Sema & Vadillo (1985). These developments raised many interesting and 
important issues of concern that are useful to practitioners in computational sciences. Some of 
the issues are: 

(a) Can recent advances in dynamical systems provide new insights into better understanding 
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of numerical algorithms and the construction of new ones? 

(b) Can these advances aid in the determination of a more reliable 'riterion on the use of 
existing numerical schemes for strongly nonlinear problems? 

(c) What is the influence of finite time steps and grid spacings rather an time steps and 
grid spacings approaching zero on the overall nonlinear behavior and stabili . af the scheme in 
terms of allowable initial data and discretized parameters? 

Since the early 1990’s, the use of dynamics to address long time behavior c numerical 
schemes for IVPs and IB VPs began to flourish. The more recent work includes the (inference 
on Dynamics of Numerics and Numerics of Dynamics (University of Bristol, dy 31 - 
August 2, 1990), the Chaotic Numerics Workshop (Deakin University, Geelong, t stralia, 
July 12-16, 1993), the Conference on Dynamical Numerical Analysis (Georgia Inst, *e of 
Technology, Atlanta, Georgia, December 14-16, 1995), and the “Innovative Time Integ. ors 
Workshop” (Center for Mathematics and computer Science, Amsterdam, November 6-8, 1 \ 

the Netherlands). These conferences were devoted almost entirely to dynamical numerk 
analysis. See the proceedings and references cited therein. See also Stuart (1994, 1995, 
Humphries (1992), Hairer et al. (1989), Aves et al. (1995), Corless (1994a,b), Dieci & Estep 
(1991), and Poliashenko & Aidun (1995). The majority of the later developments concentrated 
on long time behavior of ODE solvers using variable step size based on local error controls 
(Butcher 1987). This type of local error control enjoyed much success in controlling accuracy 
and stability for transient computations. 

The caveat is that regardless of whether finite difference (and finite volume) or finite 
element methods are employed, when time-marching approaches are used to obtain steady-state 
numerical solutions, local error controls similar to that used in ODE solvers that were designed 
for accuracy purposes are neither practical nor appropriate to use, since such local step size error 
control methods might prevent the solution from reaching the correct steady-state solutions 
within a reasonable number of iterations. We note that the standard practice of using “local 
time step’ ’ (varied from grid point to grid point with the same CFL) in time-marching to the 
steady state is not the same as the variable step size based on local error controls. 

The authors believe that the understanding of the dynamics of numerics for constant step 
size is necessary from that aspect. Besides, the study of the dynamics of ODE solvers using 
variable step size based on local error control requires a knowledge of the constant step size 
case (Aves et al. 1995). In a series of papers, Yee et al. (1991), Yee & Sweby (1994, 
1995a,b), Sweby et al. (1990, 1995), Sweby & Yee (1992, 1994), and Lafon & Yee (1991, 
1992) studied the dynamics of finite discretization for constant time steps. The examples used 
in these papers were deliberately kept simple to permit explicit analysis. The approach was to 
take nonlinear model ODEs and PDEs with known explicit solutions (the most straight forward 
way of being sure what is ‘really’ happening), discretize them according to various standard 
numerical methods, and apply techniques from discrete dynamics to analyze the behavior of the 
discretized counterparts. To set the stage for later discussion, the next few subsections discuss 



8 


the connection of dynamical systems with CFD. These subsections list some of the outstanding 
issues of numerical uncertainties in CFD in which the tools of dynamical systems theory can 
play. 


2.1. Fluid Dynamics Equations as Dynamical Systems 

Most of the available fluid dynamics and CFD related texts and reference books describe the 
Euler and Navier-Stokes equations in differential form as coupled systems of nonlinear PDEs. 
These equations are rarely classified as dynamical systems. However, fluid dynamicists are 
often interested in how the flow behaves as a function of one or more physical parameters. Of 
particular interest to fluid dynamicists is locating the critical value of the physical parameter 
where the fluid undergoes drastic, changes, in the flow behavior. Some examples are the 
prediction of transition to turbulence or laminar instability as a function of the Reynolds 
number, flow separation and stall as a function of Reynolds number and angle of attack, 
rotorcraft vibration as a function of rotation speed and flight speed, the occurrence of shock 
waves as a function of the body shape and/or Mach number, and the formation of vortices, 
flutter, and other flow phenomena as a function of the angle of attack or other physical 
parameters. One also can recast the study of admissible shock wave solutions of hyperbolic 
conservation laws as the study of the dynamics of heteroclinic orbits of a system of nonlinear 
ODEs (Shearer et al. 1987). Another application is in the area of aiding the understanding of the 
topology of flow patterns (flow visualizations) of laboratory experiments, observable physical 
phenomena and numerical data. An additional important topic for CFD is the control and 
optimization of dynamical systems. This involves the application of optimization and control 
theory to dynamical systems. Researchers are beginning to use these interdisciplinary ideas to 
study, for example, the control of turbulence, the control of vortex generation and/or shock 
waves, the control of vibration in rotorcraft, and the control of aerodynamic noise such as sonic 
boom and jet noise. 

The application of dynamical system theory to the study of spatio-temporal instabilities of 
aerodynamic and hydrodynamic flows and chaotic systems in fluid dynamics was discussed 
respectively in the 1994 and 1996 von Karman Institute for Fluid Dynamics Lecture Series. 
How the solution behaves as one or more of the system parameters is varied is precisely the 
definition of dynamical systems and bifurcation theory. According to Ian Stewart (1990) 

‘ ‘Bifurcation theory is a method for finding interesting solutions to nonlinear equations by 
tracking dull ones and waiting for them to lose stability.” 

As evident from the Third and Fourth SIAM Conference on the Application of Dynamical 
Systems, May 21-24, 1995 and May 18-22, 1997, Snowbird, Utah, presentations in treating the 
various fluid flow equations as dynamical systems have pushed these topics to the forefront of 
industrial and applied mathematical research. 
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2.2. Discrete Dynamical Systems and CFD 

When we try to use numerical methods to gain insight into the fluid physics, there is an 
added new dimension to the overall problem. Even though we freeze the physical parameters of 
the governing equations, the resulting discretized counterparts are not just a nonlinear system 
of difference equations, but are also a nonlinear but discrete dynamical system on their own. 
Depending on the scheme, the discretized counterparts usually preserve the steady states of 
the continuum. In addition, the discretized counterparts possess their own dynamics which is 
usually richer than the continuum (Mitchell & Griffiths 1985, Iserles 1988, Yee et al. 1991). 
These resulting discrete dynamical systems are a function of all of the discretized parameters 
which are not present in the governing equations. This is one of the key factors in influencing 
the numerical solution to depart from the physical one if the governing equations are strongly 
nonlinear and stiff. See Section III for an introduction. 

Of course, before analyzing the dynamics of numerics, it is necessary to analyze (or 
understand) as much as possible the dynamical behavior of the governing equations and/or the 
physical problems using theories of DEs (ODEs and PDEs), dynamical systems of DEs, and 
also physical guidelines. In fact, a knowledge of the theories of DEs, dynamical behavior of 
nonlinear DEs, and the dynamical behavior of nonlinear discrete dynamical systems (difference 
equations) is a prerequisite to the study of dynamics of numerics. In an idealized situation, if 
one knows the dynamical behavior of the governing equations, one can then construct suitable 
numerical methods for that class of dynamical systems. Consequently, spurious dynamics due 
to numerics can be minimized and that computation and analysis kept pretty much in tune. 
However, as applied scientists want to push the envelope of understanding of realistic flows and 
configurations further, dependence on the numerics takes over even though rigorous analysis 
lags behind. Starting in the late 1970’s the advances in computer power resulted in attempts 
to use CFD to replace wind tunnel experiments and use numerics to understand dynamical 
systems. The gap between computation and analysis increased. The nonlinear behavior of 
commonly used algorithms in CFD was not well understood, but at the same time applied CFD 
increased the intensity of using these algorithms to solve more complex practical problems 
where the flow physics and configurations under consideration were not understood, and were 
either too costly for or not amenable to laboratory experiments. CFD was and remains in a 
stage where computation is ahead of analysis. In other words, we usually do not know enough 
about the solution behavior of the underlying DEs in practice, and we are at the stage where the 
understanding of the dynamics of the DEs and the understanding of the dynamics of numerics 
are in tandem, and they both are rapidly growing research areas. To aid the understanding 
of the scope of the situation, first, it is important to identify all the sources of nonlinearities. 
Second, it is necessary to isolate the elements and issues of numerical uncertainties due to these 
nonlinearities. 

Sources of Nonlinearities : The sources of nonlinearities that are well known in CFD 
are due to the physics. Examples of nonlinearities due to the physics are convection, 
diffusion, forcing, turbulence source terms, reacting flows, combustion related problems, or any 
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combination of the above. The less familiar sources of nonlinearities are due to the numerics. 
There are generally three major sources: 

(a) Nonlinearities due to time discretizations — the discretized counterpart is nonlinear in 
the time step. Examples of this type are Runge-Kutta methods. It is noted that linear multistep 
methods (LMMs) are linear in the time step. See Lambert (1973) for the forms of these methods. 

(b) Nonlinearities due to spatial discretizations — in this case, the discretized counterpart can 
be nonlinear in the grid spacing and/or the scheme. Examples of nonlinear schemes are the total 
variation diminishing (TVD) and essentially nonoscillatory (ENO) schemes. See Yee (1989) 
and references cited therein for the forms of these schemes. 

(c) Nonlinearities due to complex geometries, boundary interfaces, grid generation, grid 
refinements and grid adaptations — each of these procedures can introduce nonlinearities. 

The behavior of the above nonlinearities due to the numerics are not well understood. Only 
some preliminary development is recently beginning to emerge. 


2.3. Dynamics of Numerical Approximations of ODEs vs. Time-Dependent PDEs 

Recent analyses and studies have shown that spurious numerical solutions can be inde- 
pendently introduced by time and spatial discretizations. Take the case when the ODEs are 
obtained from semi-discrete approximations of PDEs. The resulting system of ODEs contains 
more parameters (due to spatial discretizations) than those in the physical problems governed 
by ODEs. The parameters due to spatial discretizations for the semi-discrete approximation 
becomes the system parameter (instead of the discretized parameter) of the resulting system 
of ODEs. Depending on the differencing scheme, the resulting discretized counterparts of a 
PDE can be nonlinear in At, the grid spacing Ax and the numerical dissipation parameters, 
even though the PDEs have only one parameter or none. One major consideration is that one 
might be able to choose a “safe” numerical method to solve the resulting system of ODEs 
to avoid spurious stable numerical solutions due to time discretizations. However, spurious 
numerical solutions, especially spatially varying spurious steady states introduced by spatial 
discretizations in nonlinear hyperbolic and parabolic PDEs for CFD applications appear to be 
more difficult to avoid due to the use of a fixed mesh. In the case of the semi-discrete approach 
such as methods of lines or finite element methods, if spurious numerical solutions due to spatial 
discretizations exist, the resulting ODE system has already inherited this spurious feature as 
part of the exact solution of the semi-discrete case. Thus care must be taken in using the ODE 
solver computer packages for PDE applications. See Lafon & Yee (1991, 1992) and Section 
3.7 for a discussion. 


2.4. Dynamics of Time-Marching Approaches 

The use of time-marching approaches to obtain steady-state numerical solutions has been 
considered the method of choice in CFD for nearly two decades since the pioneering work of 
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Crocco (1965) and Moretti & Abbett (1966). Moretti and Abbett used this approach to solve 
the inviscid supersonic flow over a blunt body without resorting to solving the steady form of 
PDEs of the mixed type. The introduction of efficient CFD algorithms of MacCormack (1969), 
Beam & Warming (1978), Briley & McDonald (1977) and Steger (1978) marked the beginning 
of numerical simulations of 2-D and 3-D Navier-Stokes equations for complex configurations. 
It enjoyed much success in computing a variety of weakly and moderately nonlinear fluid flow 
problems. For strongly nonlinear problems, the situation is more complicated. In addition 
to the understanding of the sources of nonlinearities, it necessary to isolate all elements and 
issues of numerical uncertainties due to these nonlinearities in time-marching to the steady 
state. The following isolates some of the key elements and issues of numerical uncertainties in 
time-marching to the steady state. 

Solving an IBVP with Unknown Initial-Data : ■■ When time-marching approaches are em- 
ployed to obtain steady-state numerical solutions, a BVP is transformed into an IBVP with 
unknown initial data. The time differencing in this case acts as a pseudo time. Linearized 
stability analysis indicates that a subset of the numerical solutions for certain ranges of the 
discretized parameters and numerical boundary conditions mimic the true solution behavior 
of the governing equation. However, it is less known that there exist asymptotic numerical 
solutions (including spurious steady states) that are not solutions of the continuum inside as 
well as outside the safe regions (Yee et al. 1991, Yee & Sweby 1994, 1995a,b), depending 
on the initial data. Unlike nonlinear problems, the numerical solutions of linear or nearly 
linear problems are “independent” of the discretized parameters and initial data as long as 
the discretized parameters are inside the stability limit (or the Courant-Friedrich-Lewy (CFL) 
condition). That is, the topological shapes of these solutions remain the same within the 
stability limit and accuracy of the scheme for linear behavior. Section 3.4 illustrates the strong 
dependence of the numerical solution on initial data for nonlinear DEs. It turns out that if 
constant step sizes are used, stability, convergence rate and occurrence of spurious numerical 
solutions are intimately related to the choice of initial data (or start up solution). See Section 
3.4 for an introduction. 

Reliability of Residual Test: The deficiency of the use of residual tests in detecting the 
convergence rate and the convergence to the correct steady-state numerical solutions is now 
briefly discussed. Consider a quasilinear PDE of the form 

tit = tx*, ot, c), (2.1) 

where G is nonlinear in u, u« and u„. The values a and e are system parameters. For 
simplicity, consider a two time level and a (p + q + 1) point grid stencil numerical scheme of 
the form 


u 


n+l _ 


= «? - Ax) 


( 2 . 2 ) 


for the PDE (2.1). Note that the discussion need not be restricted to explicit methods or two 
time level schemes. Let U *, a vector representing (u^ +9 ,...,tt^,...,u^_ p ), be a steady-state 
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numerical solution of (2.2). When a time-marching approach such as (2.2) is used to solve 
the steady-state equation G(u,u.,u„,a,e) = 0, the iteration typically is stopped when the 
residual H and/or some l 2 norm of the dependent variable u between two successive iterates is 
less than a pre-selected level. 

Aside from the various standard numerical errors such as truncation error, machine round-off 
error, etc., there is a more fundamental question of the validity of the residual test and/or t 2 norm 
test. If the spatial discretization happens to produce spurious steady-state numerical solutions, 
these spurious solutions would still satisfy the residual and l 2 norm tests in a deceptively smooth 
manner. Moreover, depending on the combination of time as well as spatial discretizations, it 
is not easy to check whether (7(u*, a, e) — ► 0 even though H(U*,a,e, At, A*) — » 0, 

since spurious steady states (and asymptotes) can be independently introduced by spatial and 
time discretizations. This is contrary to the ODE case, where if u* is a spurious steady state 
numerical solution of the underlying ODE du/dt = 5(u), then S(t»") / 0. Furthermore, if a 
steady state has been reached with a rapid convergence rate for (2.2), it does not imply that the 
steady state obtained is not spurious. 

Method s Used to Accelerate Convergence Process : Methods such as iterations and relaxation 
procedures, and/or convergence acceleration methods such as conjugate gradient methods have 
been utilized to speed up the convergence process (Saad 1994). Also techniques such 
as preconditioning (Turkel 1993) and multigrid (Wesseling 1992) combined with iteration, 
relaxation and convergence acceleration procedures are commonly used in CFD. Depending on 
the type of PDEs, proper preconditioners can be established for the PDEs or for the particular 
discretized counterparts. Multigrid methods can be applied to the steady PDEs or the time- 
dependent PDEs. In either case, a combination of these methods can still be viewed as pseudo 
time-marching methods (but not necessarily of the original PDE that was under consideration). 
However, if one is not careful, numerical solutions other than the desired one can be obtained 
in addition to spurious asymptotes due to the numerics. From here on the term ‘ ‘time-marching 
approaches” is used loosely to include all of the above. It is remarked that multigrid methods 
can be viewed as the (generalized) spatial counterpart of variable time step control in time 
discretizations. 

Methods in Solving the Nonlinear Algebraic Equations From Implicit Methods'. When 

implicit time discretizations are used, one has to deal with solving systems of nonlinear 
algebraic equations. Aside from the effect of the different methods to accelerate the conver- 
gence process discussed previously, we need to know how different the dynamical behavior is 
for the different procedures (e.g., iterative vs. non-iterative) in solving the resulting nonlinear 
difference equations. See Yee & Sweby (1994, 1995a,b) and the next section for a discussion. 

Mismatch m Implicit Schemes : It is standard practice in CFD to use a simplified implicit 
operator (or mismatched implicit operators) to reduce CPUs and to increase efficiency. These 
mismatched implicit schemes usually consist of the same explicit operator but different 
simplified implicit operators. The implicit time integrator is usually of the LMM type. One 
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popular form of the the implicit operator is the so called ‘ ‘delta formulation” (Beam & Warming 
1978). The original logic in constructing this type of scheme is that the implicit operators act as 
a relaxation mechanism. However, from a dynamical system standpoint, before a steady state 
is reached, the nonlinear difference equations representing each of these simplified implicit 
operators are different from each other. They have their own dynamics as a function of the time 
step, grid spacing and initial data. They also can exhibit different types of nonlinear behavior 
if one is solving strongly nonlinear time-dependent PDEs. Before a steady state has been 
reached, during transient states, the solution procedures take different paths to get to the steady 
state, depending on the implicit operator, the time steps, initial conditions, and grid spacings 
(even with the same explicit operator). Some combination of these choices can get trapped in 
a stable spurious numerical solution. Other combinations of these choices can by pass these 
traps. Some implicit operator perform better than others, depending on the physical problem. 
Even after a steady state has been reached and the residual error of the explicit operator is zero, 
the solution can still be spurious, since a stable spurious steady state would produce a machine 
zero residual error if the spurious behavior is due to the spatial discretization. Consequently, 
these mismatched implicit operators can have different spurious dynamics and/or different 
convergence rates for the entire solution procedure. Section 4.3 describes some examples. 

Nonlinear Schemes : It is well known that all of the TVD, total variation bounded (TVB) and 
ENO schemes (see Yee (1989) or references cited therein) are nonlinear schemes in the sense 
that the final algorithm is nonlinear even for the constant-coefficient linear PDE. These 
types of schemes are known to have a slower convergence rate than classical shock-capturing 
methods and can occasionally produce unphysical solutions for certain combinations of entropy 
satisfying parameters and flux limiters (in spite of the fact that entropy satisfying TVD, TVB 
and ENO schemes can suppress unphysical solutions). See Yee (1989) for a summary of the 
subject. The second aspect of these nonlinear schemes is that even if the numerical method is 
formally of more than first-order and if the approximation converges, the rate may still be only 
first-order behind the shock (not just around the shock). This can happen for systems where 
one characteristic may propagate part of the error at a shock into the smooth domain. Sjogreen 
(1996) illustrate this phenomena with examples. See Section 4.2 for a discussion. The third 
aspect of these higher-order nonlinear schemes is their true accuracy away from shocks. See 
Donat (1994), Casper & Carpenter (1995) and Section 5.4 for a discussion. 

Schemes That are Linear vs. Nonlinear jn At : The obvious classification of time-accurate 
schemes for time-marching approaches to the steady state are explicit, implicit, and hybrid 
explicit and implicit methods. A less commonly known classification of numerical schemes 
for time-marching approaches is the identification of schemes that are linear or nonlinear in 
the time step (At) parameter space when applied to nonlinear DEs. As mentioned before, 
all LMMs (explicit or implicit) are linear in At and all multistage Runge-Kutta methods are 
nonlinear in At. Lax-Wendroff and MacCormack type of non-separable full discretizations 
also are nonlinear in At. A desirable property for a scheme that is linear in At is that, if the 
numerical solution converges, its steady-state numerical solutions are independent of the time 
step. On the other hand, the accuracy of the steady-state numerical solutions depends on At 
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if the scheme is nonlinear in At. Certain of these types of schemes are more sensitive to At 
than others. For example, Lax-Wendroff and MacCormack methods (MacCormack 1969) are 
more sensitive than the Lerat variant (Lerat & Sides 1988). A less known property of schemes 
that are nonlinear in At is that this type of scheme has an important bearing on the existence 
of spurious steady-state numerical solutions due to time discretizations. Although schemes like 
LMMs are immune from exhibiting spurious steady-state numerical solutions, as seen in Yee 
& Sweby (1994, 1995a,b), a wealth of surprisingly nonlinear behavior of implicit LMMs that 
had not been observed before was uncovered by the dynamical approach. See the next section 
for a review. 


Adaptive Time Step Bated on Locaf Error Control . It is a standard practice in CFD to use 
“local time step” (varied from grid point to grid point using the same CFL) for nonuniform 
grids. However, except in finite element methods, an adaptive time step based on local error 
control is rarely use in CFD. An adaptive time step is built in for standard ODE solver computer 
packages (Butcher 1987). It enjoyed much success in controlling accuracy and stability for 
transient (time-accurate) computations. The issue is to what extent this adaptive local error 
control confers global properties in long time integration of time-dependent PDEs and whether 
one can construct a similar error control that has guaranteed and rapid convergence to the 
correct steady-state numerical solutions in the time-marching approaches for time-dependent 
PDEs. See Section 3.6 for a discussion. 

Nonunique Steady-State Solutions of Nonlinear DEs vs. Spurious Asymptotes: The phe- 
nomenon of generating spurious steady-state numerical solutions (or other spurious asymptotes) 
by certain numerical schemes is often confused with the nonuniqueness (or multiple steady 
states) of the governing equation. In fact, the existence of nonunique steady-state solutions of 
the continuum can complicate the numerics tremendously and is independent of the occurrence 
of spurious asymptotes of the associated scheme. But, of course, a solid background in the 
theory of nonlinear ODEs and PDEs and their dynamical behavior is a prerequisite in the 
study of the dynamics of numerical methods for nonlinear PDEs. A full understanding of 
the subject can shed some light on the controversy about the “true” existence of multiple 
steady-state solutions through numerical experiments for certain flow types of the Euler and/or 
Navier-Stokes equations. 


III . Dynamics of Numerics for Elementary Examples 


With the aid of elementary examples, this section reviews the fundamentals of spurious 
behaviors of commonly used time and spatial discretizations in CFD. These examples consist 
of nonlinear model ODEs and PDEs. The numerical schemes considered for these nonlinear 
model ODEs and PDEs were selected to illustrate the following different nonlinear behavior of 
numerical methods: 

• Occurrence of stable and unstable spurious asymptotes above the linearized stability limit 



of the scheme (for constant time steps) 

• Occurrence of stable and unstable spurious steady states below the linearized stability limit 
of the scheme (for constant time steps) 

• Stabilization of unstable steady states by implicit and semi-implicit methods 

• Interplay of initial data and time steps on the occurrence of spurious asymptotes 

• Interference with the dynamics of the underlying implicit scheme by procedures in solving 
the nonlinear algebraic equations (resulting from implicit discretization of the continuum 
equations) 

• Dynamics of the linearized implicit Euler scheme solving the time-dependent equations vs. 
Newton’s method solving the steady equation 

• Spurious dynamics independently introduced by spatial and time discretizations 

• Convergence problems and spurious behavior of high-resolution shock-capturing methods 

• Numerically induced & suppressed chaos, and numerically induced chaotic transients 

• Spurious dynamics generated by grid adaptations 

Here “spurious numerical solutions (and asymptotes)” is used to mean numerical solutions 
(asymptotes) that are solutions (asymptotes) of the discretized counterparts but are not solutions 
(asymptotes) of the underlying DEs. Asymptotic solutions here include steady-state solutions, 
periodic solutions, limit cycles, chaos and strange attractors. See Thompson & Stewart (1986) 
and Hoppensteadt ( 1 993) for the definition of chaos and strange attractors. Due to the complexity 
of the nonlinear analysis, a logical breakdown of the key nonlinear studies can consist of 

• IVPs of explicit and implicit temporal discretizations 

• BVPs of linear and nonlinear spatial discretizations 

• IB VPs of time-accurate schemes 

• IB VPs of time-marching to the steady state 

• Nonlinearities introduced by grid generation, grid adaptation and complex geometries 

Except for Section 3.6, details of these examples can be found in our earlier papers. Sections 
3.1 - 3.6 discuss IVPs of temporal discretizations. Section 3.7 discusses IB VPs of temporal and 
spatial discretizations. Section 4.2.3 discusses the dynamics of grid adaptation. 

3.1. Preliminaries 

Consider an autonomous nonlinear ODE of the form 

^ = o5(«), (3.1) 

where a is a parameter and S(u) is nonlinear in u. Autonomous here means that t does 
not appear explicitly in the function 5. (For simplicity of discussion, we will not consider 
nonautonomous ODEs and autonomous ODEs that are nonlinear in a.) 

A fixed point u* of an autonomous system (3.1) is a constant solution of (3.1) or 
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S(u*) = 0. (3.2) 

The terms “equilibrium points”, “critical points”, “singular points”, “stationary points”, 
‘ ‘asymptotic solutions’ ’ (we are excluding periodic solutions for the current definition), * ‘steady- 
state solutions” and “fixed points” are sometimes used with slightly different meanings in 
the literature, for example, in bifurcation theory. However, for the current discussion and for 
the majority of the nonlinear dynamics literature, these terms are used interchangeably. Note 
that certain researchers reserve the term ‘ ‘fixed point’ ’ for difference equations (discrete maps) 
only. 

Consider a nonlinear discrete map from the finite discretization of (3.1) 

u n+1 = u n + !>(« n ,r), (3.3) 

where r = a At and D( u n , r) is linear or nonlinear in r depending on the numerical scheme. A 
similar analysis applies if (3.3) involves more than two time levels. Examples to illustrate the 
dependence on the numerical schemes for cases where D is linear or nonlinear in the parameter 
space will be given in a subsequent section. 

A fixed point u* of (3.3) (or fixed point of period 1) is defined by u n+1 = u n , or 

u* = u* + I?(u*,r); (3.4a) 

i.e., 

i)(u*,r) = 0. (3.4b) 

One can also define a fixed point of period p (or periodic solution of period p), where p is a 
positive integer by requiring that u n+i> = u n or 

u* = E p (u*,r), but u* ^ E h (u*,r) for 0 < k < p. (3.5) 

Here, E p (u*,r) means that we apply the difference operator E p times, where E(u n ,r) = 
u n -(- D{u n , r). For example, a fixed point of period 2 means u n+3 = u n or 

u' = E{E{u\r)). (3.6) 

In the context of discrete systems, the term “fixed point” without indicating the period 
means “fixed point of period 1” or the steady-state solution of (3.3). Interchangeably, we also 
use the term asymptote to mean a fixed point of any period. 

In order to illustrate the basic idea, the simplest form of the Ricatti ODE, i.e., the logistic 
ODE with 


du 

dt 


= aS(u) = ou(l — u) 


(3.7a) 
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is considered. For this ODE, the exact solution is 


u(t) = 


ti 


0 


u° + (1 — u°)e _at ’ 


(3.7b) 


where «° is the initial condition. The fixed points of the logistic equation are roots of 
u*(l — u*) = 0; it has two fixed points u* = 1 and u* = 0. 

To study the stability of these fixed points, we perturb the fixed point with a disturbance £, 
and obtain the perturbed equation 


Next, 5(u* + £) can be expanded in a Taylor series around u*, so that 


(3.8) 


ti. 

dt 


= a 


S(u*) + S u (u*)£ + -Suu(u*)£ H 


(3.9) 


where S u (u*) = ~| .. Stability can be detected by examining a small neighborhood of the 
fixed point provided that, for a given a, u* is a hyperbolic point (Seydel 1988) (i.e., if the 
real part of aS„(u*) ^ 0). Under this condition £ can be assumed small, its successive powers 
£ 2 ,£ s , ... can normally be neglected and the following linear perturbed equation is obtained 


§ = aS u («')f. (310) 

The fixed point u* is asymptotically stable if aS u (u*) < 0, whereas u* is unstable if 
a S u ( u *) > o. This stability analysis will hereafter be referred to as the linearized stability 
analysis around the fixed point tt*. If cxS u (u * ) = 0, a higher order perturbation is necessary. If 
ti* is not a hyperbolic point, the behavior of (3.10) does not infer the behavior of the original 
unperturbed equation. 

For the logistic ODE, the fixed points are hyperbolic. Thus the linearized analysis suffices 
(i.e., the original equation has the same local behavior as the perturbed equation (3.10)). If we 
perturb the logistic equation around the fixed point with a > 0, we find that u* = 1 is stable 
and u* = 0 is unstable. It is well known that the global asymptotic solution behavior of the 
logistic ODE is that for any u° > 0, the solution will eventually tend to u* = 1. Figure 3.1 
shows the asymptotic solution behavior of the logistic ODE. 

Now, let us look at three of the well known schemes for IVPs of ODEs. These are 
explicit Euler (Euler or forward Euler), leapfrog and Adams-Bashforth. For the ODE (3.1), the 
dynamical behavior of their corresponding discrete maps is well established. The explicit Euler 
method is given by 


„"+! =t," + r5(u n ), 


(3.11) 
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and after a linear transformation it is the well known logistic map (Hoppensteadt 1993). The 
leapfrog method can be written as 

u n+1 = u”" 1 + 2r5(«"), (3.12) 

and it is a form of the Henon map (Devaney 1987). The Adams-Bashforth method yields 

u n+i = u n + ' [ 3 5(« n ) - 5(u n-1 )l , (3.13) 

again a variant of the Henon map that has been discussed by Proffer (1985) in detail. 

We can determine fixed points of the discrete maps (3. 1 1 )-(3. 13) and their stability properties 
in a manner similar to that for the ODE. It turns out that all three of the discrete maps have the 
same fixed points as the ODE (3.1) — a desired property which is important for obtaining the 
correct steady states of nonlinear DEs numerically. An examination of (3.1 1)-(3.13) reveals 
that the discretized parameter r appears linearly in these discrete maps. If nonlinearity in the 
parameter space At is introduced into the discretized counterpart, it increases the possibility 
that the dynamics of numerics deviates from the dynamics of the continuum. As can be seen 
in Yee et al. (1991), a necessary condition for the occurrence of spurious steady states by any 
scheme is that the discretized parameter should appear nonlinearly in the underlying discrete 
map. Consequently the existence of spurious steady states is not possible for (3.1 1)-(3.13). 

The corresponding linear perturbed equation for the discrete map (3.3), found by substituting 
u n = u* + in (3.3) and ignoring terms higher than £ n , is 

( n+1 = t n [1 + AtD u (u\ Af)] • (3-14) 

Here the parameter a of the ODE has been absorbed in the parameter At based on the 
assumption that a does not appear explicitly in S(u). Depending on the scheme, D( u n , At) 
might be nonlinear in At. 

For stability we require 

|1 + A<D u (u*,Ai)| < 1. (3.15) 

Again, for |1 + AtD u (u*, At)| = 1, a higher order perturbation is necessary. For a fixed point 
of period p the corresponding linear perturbed equation and stability criterion are 

r +p = (" E p (u .At), (3.16a) 

and 


with 


|^(u%A<)|<1, 


(3.16b) 
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E£(u n ,A<) = ^-E(u n+p -\At)...-^E(u n ,At). (3.16c) 

For S(u) = it(l — u), the stability of the stable fixed points of period 1 and 2 for discrete 
maps (3. 1 1 )-(3. 13) with r = a At are 

Explicit Euler: 

u* = 1 stable if 0 < r < 2 

period 2 stable if 2 < r < \/6 

Leapfrog: 

u* — 1 unstable for all r-> 0 

chaotic solutions exist for all r no matter how small 

Adams-Bashforth : 

u* = 1 stable if 0 < r < 1 

period 2 stable if 1 < r < V2. 

Figure 3.2a shows the stable fixed point diagram of period 1,2,4, 8 obtained for the 
explicit Euler scheme by solving numerically the roots of (3.11) (by setting u n+1 = u n ) for 
S(u) = u(l — u). The r axis is divided into 1000 equal intervals. The numeric labeling of the 
branches denotes their period. The subscript 4 4 E’ ’ on the period 1 branch indicates the stable 
fixed point of the DE. From here on, r = aAt where a = a in all of the plots. The change of 
notation inside the plots is due to the plotting package. 

All of these three examples share a common property of not exhibiting spurious steady 
states. Two of these three examples serve to illustrate that operating with a time step beyond the 
linearized stability limit of the stable fixed points of the nonlinear ODEs does not always result 
in a divergent solution; spurious asymptotes of higher period can occur. This is in contrast to 
the ODE solution, where only a single stable asymptotic value u* = 1 exists for any a > 0 and 
any initial data u° > 0. The spurious asymptotes, regardless of the period, stable or unstable, 
are solutions in their own right of the discrete maps resulting from a finite discretization of the 
ODE. 


3.2. Spurious Asymptotic Numerical Solutions for Constant Time Steps 


For the previous section we purposely picked the type of schemes that do not exhibit spurious 
fixed points, but allow spurious fixed points of period higher than 1 . In this section numerical 
methods are purposely chosen so that the discretized parameter appears nonlinearly in the 
underlying discrete maps. Consequently, existence of spurious steady-state numerical solutions 
in these examples is possible. 
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3.2.1. Explicit Time Discretizations 

Consider two second- and third-order Runge-Kutta method, namely, the modified Euler 
(R-K 2) and the improved Euler (R-K 2), Heun (R-K 3), Kutta (R-K 3); also consider the 
fourth-order Runge-Kutta method (R-K 4), and two predictor-corrector methods (P-C 2 and 
P-C 3) (Lambert 1973) of the forms 

Modified Euler (R-K 2) method: 

u n+1 =u n + rS (u n + ^5 n ), 5" = S(u"), (3.17) 

Improved Euler (R-K 2) method: 

ti n+I =u n + ^ S n + S(u n + rS n ) , (3.18) 


Heun (R-K 3) method: 



Kutta (R-K 3) method: 



(3.19) 


(3.20) 
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R-K 4 method: 


u 


"+1 _ 


« + 


— + 2fej + k^j 


k\ =5 
k 2 


=*(“"+K) 

*> = s (“”+i‘ i ) 

- = 5" (u n + , 


Predictor-corrector for m = 2,3 (P-C 2 & P-C 3): 


(3.21) 


«(°) = u" + rS r 

u (fc+i) = u n + r 

2 

«- +1 = « n + ? 


5 n + 5 (fc) 

5 n + 5 (m_1) 


It = 0, — 1 


(3.22) 


3.2.2. Fixed Point Diagrams 


Using the same procedures as before, one can obtain the fixed points for each of the above 
schemes (3.17) - (3.22). Figures 3.2b - 3.2f show the stable fixed point diagrams of period 
1,2,4 and 8 for selected schemes for S(ti) = tx(l — u). The unstable fixed points of any period 
are not plotted. See Yee et al. (1991) for the unstable fixed point diagrams. Some of the fixed 
points of lower period were obtained by closed form analytic solution and/or by a symbolic 
manipulator such as MAPLE (1988) to check against the computed fixed point. The majority 
were computed numerically. The stability of these fixed points was examined by checking the 
discretized form of the appropriate stability conditions. The domain is chosen so that it covers 
the most interesting part of the scheme and ODE combinations, and is divided into 1 000 equal 
intervals. In other words, spurious asymptotes may occur in other parts of the domain as well. 
The numeric labeling of the branches denotes their period, although some labels for period 4 
and 8 are omitted due to the size of the labeling areas. Again, the subscript E on the main 
period one branch indicates the stable fixed point of the ODE while the subscript “S” indicates 
the spurious stable fixed points introduced by the numerical scheme. Spurious fixed points of 
period higher than one are obvious (since the ODEs under discussion only possess steady-state 
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solutions) and are not labeled with a subscript “S”. Note that these diagrams, which for the 
most part appear to consist of solid lines, actually consist of points, which are only apparent 
in areas with high gradients. 

To contrast the results, similar stable fixed point diagrams were also computed for the ODE 

— = au(l — u)(5 — u), 0 < 6 < 1, (3.23) 

dt 

that is, for a cubic nonlinearity for S(u ) = u(l — u)(6 — u). The stable fixed point for the ODE 
(3.23) in this case is u* = b and the unstable ones are u* = 0 and u* = 1. For any 0 < u° < 1 
and any a > 0, the solution will asymptotically approach the only stable asymptote of the ODE 
u* = b. 

By looking at the roots of the underlying discrete maps of (3.17)-(3.22), it is readily realized 
that r appears nonlinearly in these discrete maps. In fact, the maximum number of stable and 
unstable fixed points (real and complex) for each of the studied schemes (3.17)-(3.22) varied 
from 4 to 16 for S(u) = u(l — u) and 9 to 81 for 5(u) = u(l — u)(b — u), depending on 
the numerical method and the r value. Aside from the real steady states, these roots might be 
unstable and/or complex for certain r values but not for others. Fig. 3.3 shows the stable fixed 
point diagram by the modified Euler for four different values of b (the unstable fixed points are 
excluded from the plots). 

Aside from the striking difference in topography in the stable fixed point diagrams of the 
above methods and ODE combinations, all of these diagrams have one common feature: they 
all exhibit spurious stable fixed points of period higher than one. In the majority of cases, 
they also exhibit stable spurious steady states. In some of the instances, these spurious fixed 
points are outside the interval of the stable and unstable fixed points of the ODEs. Others not 
only lie below the linearized stability limit but also in the region between the fixed points of 
the ODEs and so could be very easily achieved in practice. For example, in Fig. 3.2b, the 
modified Euler scheme for the logistic ODE, the linearized stability limit of period 1^ is r = 2. 
But depending on the value of r, two stable fixed points of period 1 (one is spurious) can 
exist at the same time for 0 < r < 1.236. For the R-K 4 method applied to the logistic ODE, 
one can see from Fig. 3.2d that spurious steady states which exist for 2.75 < r < 2.785 are 
below the linearized stability limit of the 1 e branch. For the modified Euler method applied 
to du/dt = au(l — u)(5 — u), it is interesting to see the changing behavior of stable spurious 
steady states as the stable fixed point u* = b is varied between 0 and 0.5. 

A unified analysis of the above for the standard explicit Runge-Kutta methods is reported in 
Griffiths et al. (1992a). Tables 3.1 - 3.4, taken from Griffiths et al. (1992a), show the true and 
spurious asymptotes of selected schemes. Some entries are marked with an asterisk to indicate 
where stable fixed points are known to exist but no closed analytic form has yet been found. 

Historically, Iserles (1988) was the first to show that while LMMs for solving ODEs possess 
only the fixed points of the original ODEs, popular Runge-Kutta methods may exhibit spurious 
fixed points. Iserles et al. (1990) and Hairer et al. (1989) classified and gave guidelines 
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and theory on the types of Runge-Kutta methods that do not exhibit spurious period one or 
period two fixed points. Humphries (1991) showed that under appropriate assumptions if stable 
spurious fixed points exist as the time-step approaches zero, then they must either approach 
a true fixed point or become unbounded. Hence repeating the integration with a smaller step 
size will ultimately make the spurious behavior apparent. However, convergence in practical 
calculations involves a finite time step At that is not small as the number of integrations 
n — ► oo rather than At — » 0, as n — ► oo. The work in Yee et al. (1991), Yee & Sweby (1994, 
1995a,b), and Lafon & Yee (1991, 1992), Sweby et al. (1990, 1995), Sweby & Yee (1991), 
and Griffiths et al. (1992) attempted to provide some of the global asymptotic behavior of time 
discretizations when finite fixed but not extremely small At is used. Vadillo (1997) relates 
existence of spurious steady states with the numerical solution of the exact steady state near 
At — ► oo. As can be seen later, stable and unstable spurious fixed points of all orders need to 
be accounted for in the study of spurious behavior of numerical schemes. 


3.3. Bifurcation Diagrams 

This section discusses another method for obtaining the stable fixed point diagrams or 
bifurcation diagrams before illustrating the symbiotic relationship between permissibility of 
spurious steady states and initial data in fixed time step computations. 

“Full” Bifurcation Diagram (" Complete Fixed Point Diagram”): If one obtains the full 
spectrum of these fixed points of any order as a function of the step size, the fixed point 
diagram is sometimes referred to as the “full” bifurcation diagram. In other words, the “full” 
bifurcation diagram exhibits the complete asymptotic solutions of the discretized counterparts 
as a function of the discretized parameter r. In computing the “full” bifurcation diagram, 
searching for the roots and testing for stability of highly complicated nonlinear algebraic 
equations (for fixed points of higher period and/or complex nonlinear DEs combination) can 
be expensive and might lead to inaccuracy. In certain instances, one might be able to obtain 
the bifurcation diagram by some type of continuation method. The most popular one is called 
the pseudo arclength continuation method and was devised by Keller (1977). However, the 
majority of the continuation type methods require known start up solutions for each of the 
main bifurcation branches before one can continue the solution along a specific main branch. 
For problems with complicated bifurcation patterns, the arclength continuation method cannot 
provide the complete bifurcation diagram without the known start up solutions. In fact, it is 
usually not easy to locate even just one solution on each of these branches, especially if spurious 
asymptotes exist. 

Computed Bifurcation Diagram: A numerical approach for obtaining a “computed” bifur- 
cation diagram (not necessarily the full bifurcation diagram, as explained later) of the resulting 
discretized counterpart consists of iteration of the underlying discrete map. In other words, 
this type of computed bifurcation diagram for the one-dimensional discrete map displays the 
iterated solution tt n vs. r after iterating the discrete map for a given number of iterations with a 
chosen initial condition for each of the r parameter values. For the figures shown later, with a 
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given interval of r and a chosen initial condition, the r axis is divided into 500 equal spaces. In 
each of the computations, the discrete maps were iterated with 600 preiterations (more or less 
depending on the DE and scheme combinations) and the next 200 iterations were over plotted 
in the same diagram for each of the 500 r values. The preiterations were necessary in order for 
the solutions to settle to their asymptotic value. A high number of iterations were overlaid on 
the same plot in order to detect periodic orbits (in this case periods of up to 200) or invariant 
sets. The reader is reminded that with this method of computing the bifurcation diagrams, only 
the stable branches are obtained. The domains of the r and u n axes are chosen to coincide with 
the stable fixed point diagrams shown previously. As explained later, even though all of these 
discrete maps possess periodic solutions of period n for arbitrarily large n and stable chaotic 
solutions, no attempt was made to compute all of the spurious orbits of any order or chaotic 
solutions. The purpose of the present discussion is to show the spurious behavior and these 
computations suffice to serve the purpose. 

Examples : Figure 3.4 shows the bifurcation diagram of the explicit Euler method applied to 
the logistic ODE with an initial condition u° = 0.5. It is of interest to know that in this case 
the bifurcation diagram looks practically the same for any u° > 0. This is due to the fact that 
no spurious fixed points exist for r < 2 and no spurious asymptotes of low period exist for 
r < 2.627. One quickly observes that using the arclength continuation method for this discrete 
map is the most efficient way to obtain its bifurcation diagram. However, this is not the case for 
other methods to be discussed later. Comparing the bifurcation diagram with Fig. 3.2a, one can 
see that if we had computed all of the fixed points of period up to 200 for Fig. 3.2a, the resulting 
fixed point diagram would look the same as the corresponding bifurcation diagram (assuming 
600 iterations of the logistic map are sufficient to obtain the converged stable asymptotes of 
period up to 200 and the chosen initial data are appropriate to cover the basins of all of the 
periods in question). The numeric labeling of the branches in the bifurcation diagram denote 
their period, with only the essential ones labeled for identification purposes. 

The noise appearing on the 1 e branch near the bifurcation point r = 2 of the linearized 
stability limit of the fixed point u* = 1 indicates that 600 iterations of the logistic map are not 
sufficient to obtained the converged stable asymptotes. This phenomenon is common to other 
bifurcation points of higher periods as well as the rest of the corresponding bifurcation for the 
studied schemes. See Yee et al. (1991), and Yee & Sweby (1994, 1995a, b) for additional 
details. In fact, the slow convergence of using a time step that is near the linearized stability of 
the scheme (bifurcation point) might be due to this fact. 

We note that the explicit Euler applied to the logistic ODE resulted in the famous logistic 
map. Unlike the underlying logistic ODE, it is well known that the logistic map possesses very 
rich dynamical behavior such as period-doubling (of period 2” for any positive n) cascades 
resulting in chaos (Feigenbaum, 1978). One can find Fig. 3.4 appearing in most of the 
elementary dynamical systems text books. The exact values of r for all of the period-doubling 
bifurcation points and chaotic windows (intervals of r) were discovered by Feigenbaum in the 
late 1970’s. Interested readers should consult these elementary text books for details. In other 
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words, one can obtain the analytical (exact) behavior of the spurious asymptotes and numerical 
(spurious) chaos of the logistic map. The next section explains why using a single initial datum 
in computing the bifurcation diagrams for schemes that exhibit spurious asymptotes does not 
necessarily coincide with the fixed point diagram (or full bifurcation diagram). It is interesting 
to note that the corresponding bifurcation diagrams of the respective discrete maps produced 
by the remaining studied schemes consist of unions of “logistic-map-like” bifurcations and/or 
‘ ‘inverted logistic-map-like’ ’ bifurcations with similar yet slightly complicated period-doubling 
cascades resulting in chaos. See Section 3.4 for additional discussions. 

Types of Bifurcations: In all of the fixed point diagrams shown previously, the majority of 
the bifurcation phenomena can be divided into three kinds; these are flip, supercritical, and 
transcritical bifurcations (Seydel 1988). Figure 3.5a shows the schematic of typical types of 
steady bifurcations (bifurcations of' asymptotes other than steady ones are not shown). Figure 
3.5b shows examples of these three types of bifurcation for the logistic ODE using the modified 
Euler, improved Euler and R-K 4 methods. Figure 3.5b also shows a comparison of the stable 
and unstable fixed points of periods 1 and 2. Although the modified Euler and the R-K 4 
methods experience a transcritical bifurcation, they have different characteristics. See Fig. 3.5a 
for the different types of transcritical bifurcations. Note that the flip bifurcation looks very 
much like the supercritical (steady) bifurcation. However, in the flip bifurcation, the solution 
becomes periodic after the flip bifurcation point. 

For the bifurcation of the first kind, the paths (spurious or otherwise) resemble period 
doubling bifurcations (flip bifurcation) similar to the logistic map. See Figs. 3.2a,e (for r — 2) 
for examples. The second kind is a steady or supercritical bifurcation. It occurs most often at the 
main branch 1 e with the spurious paths branching from the correct fixed point as it reaches the 
linearized stability limit, quite often even bifurcating more than once as r increases still further 
before the onset of period doubling bifurcations. See Figs.3.2c,f (for r = 2) for examples. 
Using the P-C 3 method to solve (3.7) , more than one consecutive steady bifurcation occurs 
before period doubling bifurcations. Follow the Is labels on Fig. 3.2f. Although figures are not 
shown for ODE (3.23) with b = 0.5, the improved Euler experiences two consecutive steady 
bifurcations before period doubling bifurcation occurs (see Yee et al. 1991, for details). Using 
the P-C 3 method to solve (3.23), four consecutive steady bifurcations occur before period 
doubling bifurcations. The modified Euler and R-K 4 methods, however, experience only one 
steady bifurcation before period doubling bifurcations occur. 

The third kind of bifurcation again occurs most often at the main branch 1 e - The spurious 
paths near the linearized stability limit of 1® experience a transcritical bifurcation. This is 
another kind of steady bifurcation. See Fig. 3.2b (for r = 2), Fig. 3.2d (for r near 2.75), 
Fig. 3.2f (for r near 3.4) and Fig. 3.3 (follow the 1® branch) for examples. Notice that the 
occurrence of transcritical and supercritical bifurcations is not limited to the main branch \ e - 
See Fig. 3.3 for examples. At the stability limit of the true fixed point, only the modified Euler 
and R-K 4 undergo transcritical bifurcation. 

As can be seen, the occurrence of flip and supercritical bifurcations is more common. In fact. 
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most of the bifurcation points shown in previous figures are of these types. The other commonly 
occurring bifurcation phenomenon is the subcritical bifurcation which was not observed in our 
two chosen 5(u) functions. With a slight change in the form of our cubic function S(tx), a 
subcritical bifurcation can be achieved (Seydel 1988). See elementary text books on bifurcations 
of discrete maps (Seydel 1988) for a discussion of these four types of bifurcation phenomena. 
A consequence of the latter three bifurcation behaviors is that bifurcation diagrams computed 
from a single initial condition u° will appear to have missing sections of spurious branches, 
or even seem to jump between branches. This is due to the existence of spurious asymptotes 
of some period or more than one period, and its dependence on the initial data. Section 3.4 
discusses this issue in more detail. First, we would like to look at convergence rates that are 
near a bifurcation point. 

Slow Convergence Near Bifurcation Points : As discussed previously, the number of itera- 
tions for the computed asymptotes that are near or at the bifurcation point can be orders of 
magnitude higher than away from the bifurcation points. In fact, depending on the type of 
bifurcation and initial data, one might experience slow convergence using a time step that is 
near the linearized stability of the scheme (bifurcation points of the above four types). See Yee 
& Sweby (1994, 1995a,b) for some examples. In the worst case scenario, if the bifurcation is of 
the transcritical or subcritical type and the time step is within that range, the numerical solution 
can get trapped in a spurious steady state or a spurious limit cycle, causing nonconvergence of 
the numerical solution. 


3.4. Strong Dependence of Solutions on Initial Data (Numerical Basins of Attraction) 

Computing Bifurcation Diagrams Using A Single Initial Datum : Figures 3.6a - 3.6c show 
the bifurcation diagram of the modified Euler method for the logistic ODE with three different 
starting initial conditions (I.C.). In contrast to the explicit Euler method as shown in Fig. 3.4, 
none of these diagrams look alike. One can see the influence and the strong dependence of 
the asymptotic solutions on the initial data. For certain initial data and At value combinations, 
spurious dynamics can be avoided. Yet for other combinations, one can never get to the correct 
steady state. In other words, it is possible that for the same At but two different initial data 
or vice versa, the scheme can converge to two different distinct numerical solutions of which 
one or neither of them is the true solution of the underlying ODE. Thus in a situation where 
there is no prior information about the exact steady-state solution, and where a time-marching 
approach is used to obtain the steady-state numerical solution when initial data are not known, 
a stable spurious steady-state could be computed and mistaken for the correct steady-state 
solution. Figure 3.6d shows the corresponding “full” bifurcation diagram, their earlier stages 
resembling the fixed point diagram 3.2b which showed only fixed points up to period 8. See 
Yee et al. (1991) for an example where overplotting a number of initial data, but not the 
appropriate ones, is not sufficient to cover all of the essential spurious branches. The strong 
dependence of solutions on initial data is evident from the various examples in which this type 
of behavior is present. We note that if one uses the pseudo arclength continuation type method 
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without solving for the roots of the spurious fixed point, one only knows one starting solution 
(i.e., the exact steady states of these two ODEs). The continuation method, in this case, only 
produces the branch of the bifurcation diagram originating from the 1^ branch of the curve. 

Computed Full Bifurcation Diagram : In order to compute a ‘ ‘full’ ’ bifurcation diagram using 
this numerical approach, we must overplot all of the individual bifurcation diagrams of existing 
asymptotes of any period and chaotic attractors obtained by using the entire domain of u values 
as starting initial data. Thus, a better method in numerically approximating the full bifurcation 
diagram is dividing the domain of interest of the u axis into equal increments and using these 
v values as initial data. The “full” bifurcation diagram is obtained by simply overplotting all 
of these individual diagrams on one. Figs. 3.7 and 3.8 show the “full” bifurcation diagrams 
for the corresponding fixed point diagrams shown previously. Note that the full bifurcation 
diagram computed this way might miss-some of the windows of- bifurcations that occur inside 
the intervals of the adjacent r and/or the initial data values. 

It is noted that for the cases when one knows the bifurcation pattern of a specific discrete map, 
the actual number of the initial data points used for that full bifurcation diagram computations 
do not have to completely cover the entire domain of u as long as these initial data cover all of 
the basins of attraction of the asymptotes (i.e., which initial data lead to which asymptotes). See 
the next section for a definition and discussion. That is, at least one initial data point is used 
from each of the basins of the asymptotes. No attempt has been made to compute the complete 
full bifurcation diagram, since this is very costly and involves a complete picture of the existing 
asymptotes of any period and chaotic attractors for the domain of interest in question. See 
remarks in Section 3.2 on computed bifurcation diagrams for an explanation. Here, we use the 
term “full bifurcation diagram” to mean “computed bifurcation diagram with sufficient initial 
data to cover the essential lower order periods”. Without loss of generality, from here on we 
use the word bifurcation diagram to mean the computed (and approximated) full bifurcation 
diagram. 

From Figs. 3.7 and 3.8, one can conclude that all of the studied explicit methods eventually 
undergo “period doubling bifurcations” leading to the “logistic-map-type bifurcations”. 
The term “logistic-map-type bifurcations” here means that the behavior and shape of the 
bifurcations resemble the logistic map as shown in Fig. 3.4. The ranges of the r values in 
which logistic-map-type bifurcations occur are not restricted to the 1.® branch of the bifurcation 
diagram. The birth of the logistic-map-type bifurcations can occur below or beyond the 
linearized stability limit of the true steady state of the governing equation. To aid the reader, 
Figs. 3.7 and 3.8 indicate the major stable fixed points of periods up to 4. Basically most of 
the Is branches of the bifurcations are logistic-map-type. For example, the modified Euler 
method experiences two period doubling bifurcations for the ODE (3.7). For the ODE (3.23), 
the modified Euler method experiences at least two to three period doubling bifurcations, 
depending on the b values. For other methods, the situation is slightly more complicated. 

Besides the regular logistic-map-type bifurcations, some of these methods undergo the so 
called “inverted logistic-map-type” of period doubling bifurcations. The shape of this type 
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of bifurcation resembles the reverse image of Fig. 3.4. See, for examples. Fig. 3.7b for 
1.62 < r < 1.67 and Fig. 3.7c for 3.15 < r < 3.3. 

The above example explains the role of initial data in the generation of spurious steady-state 
numerical solutions, stable and unstable spurious numerical chaos and other asymptotes. Section 
3.5 illustrates the role of initial data in the occurrence of stabilizing unstable steady states of 
the governing equation and the introduction of stable and unstable spurious numerical chaos 
and other asymptotes by implicit LMMs. Next, how basins of attraction can complement the 
bifurcation diagrams in gaining more detailed global asymptotic behavior of time discretizations 
for nonlinear DEs will be illustrated. 

“Exact” Basin of Attraction vi. “Numerical” Batin of Attraction : Associated with an 

asymptote, the basin of attraction of _an asymptote (for the DEs or their discretized coun- 
terparts) is a set of all initial data asymptotically approaching that asymptote. We use the terms 
“exact” basin of attraction and “numerical” basin of attraction to mean the basin of attraction 
of the DE and basin of attraction of the underlying discretized counterpart. Although all of 
the numerical basins of attraction shown later are obtained numerically, the term * ‘numerical 
basins of attraction” here pertains to the computed basins of attraction of all of the asymptotes 
for the underlying discrete map. 

For the logistic ODE, the “exact” basin of attraction for the only stable fixed point u* = 1 
is the entire positive plane for all values of a > 0. The basin of attraction for the ODE (3.23) 
for the stable fixed point u* = 6is0<ti<l for allo>0,0<6<l. However, the 
situation for the corresponding numerical basins of attraction for the various schemes is more 
complicated, since each asymptote, stable or unstable, spurious or otherwise, possesses its own 
numerical basin of attraction. Intuitively, one can see that in the presence of stable and unstable 
spurious asymptotes, the basins of the true stable and unstable steady states (and asymptotes 
if they exist) can be separated by the numerical basins of attraction of the stable and unstable 
spurious asymptotes. Consequently, what is part of the true basin of a particular fixed point of 
the governing equation might become part of the basin of the spurious asymptotes. For implicit 
methods that can stabilize unstable steady states of the governing equation (to be discussed 
later), the basin of attraction associated with the particular stabilizing steady state can consist 
of the union of parts of basins from other true asymptotes. In other words, the allowable 
initial data of the numerical method associated with a particular asymptote of the DE can 
experience enlargement or contraction, can become null or consist of a union of disjoint regions. 
These regions can be fractal like. Therefore, keeping track of which initial data lead to which 
asymptotes (exact or spurious) of the underlying discrete maps becomes more complicated as 
the number of spurious asymptotes increases. 

On the other hand, the computed bifurcation diagram cannot distinguish between the types 
of bifurcation and the periodicity of the spurious fixed points of any order. With the numerical 
basins of attraction and their respective bifurcation diagrams superimposed on the same plot, the 
type of bifurcation and the determination of which initial data lead to which stable asymptotes, 
become apparent. In order to obtain the corresponding numerical basins of attraction for the 
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schemes discussed above, one immediately realizes that, in most cases, a numerical approach 
is the only recourse until more theoretical tools for searching for the basin boundaries of 
general discrete maps become available. We would like to add that there are isolated theories 
or approximate methods to locate some basin boundaries for simple discrete maps or special 
classes of discrete maps. Even in this case, these methods are neither practical, nor are there 
fixed guidelines for the actual implementation of discrete maps for more complex ones of 
similar type. See Hsu (1987) for an approximate method and Friedman (1995) for numerical 
algorithms which compute connecting orbits. 

Computing Numerical Batins of Aiitaction on the Connection Machine'. The nature of 
this type of computation, especially for systems of nonlinear ODEs or PDEs, requires the 
performance of a very large number of simulations with different initial data; this can be 
achieved efficiently by the “use of the highly parallel Connection Machine (CM-2 or CM-5) 
or the IBM SP2 machine, whereby each processor could represent a single initial datum and 
thereby all the computations can be done in parallel to produce detailed global stability behavior 
and the resulting basins of attraction. With the aid of highly parallel Connection Machines, we 
were able to detect a wealth of the detailed nonlinear behavior of these schemes for systems of 
ODEs and PDEs which would have been overlooked had isolated initial data been chosen on 
the Cray or other serial or vector machine. 

Figure (3.9) shows the bifurcation diagram with numerical basins of attraction superimposed 
for the logistic ODE for four Runge-Kutta methods. (Even though one might not need to use a 
highly parallel machine to compute the basins of attraction for the scalar ODE, nevertheless this 
figure was computed on the CM-5, requiring only a few minutes for each plot). The major work 
on the CM-5 coding is on the efficient handling of data for plotting. The same plots would have 
required many orders of magnitude more CPU time on a serial or vector machine. To obtain 
a bifurcation diagram with numerical basins of attraction superimposed using the CM-5, the 
preselected domain of initial data and the preselected range of the A t parameter are divided into 
512 equal increments. For the bifurcation part of the computations, the discretized equations 
are, in general, preiterated 3000 steps for each initial datum and A t before the next 1000 
iterations (more or less depending on the problem and scheme) are plotted. The preiterations 
are necessary in order for the trajectories to settle to their asymptotic value. The high number 
of iterations are plotted (overlaid on the same plot) to detect periodic solutions. The bifurcation 
curves appear on the figures as white curves, white dots and dense white dots. While computing 
the bifurcation diagrams it is possible to overlay basins of attraction for each value of A t 
used. For the numerical basins of attraction part of the computations, with each value of A t 
used, we keep track where each initial datum asymptotically approaches, and color code each 
basin (appearing as a multi-color vertical line) according to the individual asymptotes. Black 
regions denote basins of divergent solutions. While efforts were made to match color coding 
associated with a particular asymptote of adjacent multi-color vertical lines on the bifurcation 
diagram (i.e., from one At to the next), it was not always practical or possible. Care must 
therefore be taken when interpreting these overlays. In an idealized situation it is best if we also 
know the critical value of At for the onset of unstable spurious asymptotes. However, with 
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the current method of detecting the bifurcation curve only the stable ones are detected. For 
example, a steady bifurcation would break the domain immediately after the bifurcation point 
into two different color domains, whereas the domain remains the same color immediately after 
a periodic doubling bifurcation. 

Example t : Any initial data residing in the green region in Fig. 3.9 for the modified Euler 
method belongs to the numerical basin of attraction of the spurious (stable) branch emanating 
from u = 3 and r = 1. Thus, if the initial data is inside the green region, the solution can never 
converge to the exact steady state using even a small fixed but finite At (all below the linearized 
stability limit of the scheme). Although not shown, we have computed the bifurcation diagram 
with wider ranges of r and initial data. In fact, the green region actually extends upward as 
r decreases below 1. Thus for a small range of r values very near zero, the entire domain is 
divided into two basins (not shown). As r increases, the domain divide into more than two 
basins (instead of the two for the ODE). But of course higher period spurious fixed points exist 
for other ranges of r and more basins are created within the same u domain. For r near 2 (i.e., 
near the linearized stability limit of the true steady state «* = 1, the bifurcation is transcritical. 
Using an r value slightly bigger than 2 will lead to the spurious steady state until r increases 
beyond 3.236. Consequently, there are large ranges of r below and beyond the linearized 
stability limit of the true steady state for which spurious dynamics occur. Observe the size and 
shape of these basins as r varies. The majority of these basin boundaries are fractal like. 

A similar situation exists for the R-K 4 method (Fig. 3.9), except that now the numerical 
basins of attraction of the spurious fixed points occur very near the linearized stability limit 
of the scheme, with a small portion occurring below the linearized stability limit. Although 
both the modified Euler and the R-K 4 methods experience a transcritical bifurcation for the 
logistic ODE, the transcritical bifurcation for the R-K 4 is more interesting. See Fig. 3.5 for the 
distinction between the two transcritical bifurcations. In contrast to the improved Euler method, 
the green region represents the numerical basin of attraction of the upper spurious branch of 
the transcritical bifurcation. The bifurcation curve directly below it with the corresponding red 
portion is the basin of the other spurious branch. With this way of color coding the basins of 
attraction, one can readily see (from the plots) that for the logistic ODE (3.7), the improved 
Euler method experiences one steady bifurcations before period doubling bifurcation occurs. 
The modified Euler and R-K 4 methods, however, experience a transcritical bifurcation before 
period doubling bifurcations occur. The Kutta method experiences period doubling bifurcation 
at the linearized stability limit. One way to detect the steady bifurcation from these plots is 
to look for a separate color associated with each branch of the associated bifurcation. A similar 
interpretation holds for certain types of transcritical bifurcation. See the R-K 4 plot in Fig. 3.9. 
One way to detect the period doubling bifurcation from these plots is to look for no change in 
colors associated with each branch of the associated bifurcation. For subcritical bifurcation, it 
is slightly more complicated. 

The above discussion shows the interplay between initial data, step size and permissibility of 
spurious asymptotes. It indicates that it is not just the occurrence of stable spurious numerical 
solutions that causes difficulty. Indeed such cases may be easier to detect. These spurious 
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features of the discretizations often occur but are unstable; i.e., they do not appear as an actual 
(spurious) solution because one usually cannot obtain an unstable asymptotic solution by mere 
forward time integration. However, far from being benign, they can have severe detrimental 
effects on the allowable initial data of the true solution for the particular method, hence 
causing slow convergence or possibly even nonconvergence from a given set of initial data, 
even though the data might be physically relevant. 

Due to space limitations, interested readers are referred to Yee & Sweby (1994, 1995a) for 
results for four 2x2 systems of nonlinear model ODEs. Classification of fixed points of systems 
of equations are more involved than of the scalar case. See elementary dynamical systems text 
books for details. The corresponding bifurcation diagrams and numerical basins of attraction of 
these schemes are even more involved. New phenomena exist that are absent in the scalar case. 
For example, in the. presence of multiple steady states^even for explicit methods, depending on 
the step size and initial data combination, the associated numerical basins of attraction for a true 
steady state might experience an enlargement of their basins at the expense of a contraction of 
the other asymptotes. For the scalar case, this is only possible for superstable implicit methods 
(see next section). See Yee & Sweby (1994, 1995a,b) for more details. Another example is 
that the fixed point can change type as well as stability (e.g., from a saddle point to a stable or 
unstable node). We note that for systems beyond 3 x 3, it is impossible to conduct the type of 
detailed analysis shown above. 

3.5. Global Asymptotic Behavior of Superstable Implicit LMMs 

This section reviews the superstable property of some implicit LMMs and summarizes their 
global asymptotic nonlinear behavior using the dynamical approach. Recall that the underlying 
discrete maps from using LMMs are linear in the discretized parameter At and they are exempt 
from spurious steady states. As can be seen later, the combination of implicit LMMs and the 
superstable property produce asymptotic behavior that is very different from schemes that were 
studied in the previous section. One distinct property of these types of schemes is that they 
can stabilize unstable steady states of the governing equation. Another property is that the 
numerical basins of attraction of the stable steady state can include regions of the basins of the 
unstable steady states of the governing equation. 

3.5.1. Super-stability Property 

Dahlquist et al. (1982) first defined super-stability in ODE solvers to mean the region 
of numerical stability that encloses regions of physical instability of the true solution of the 
ODE. Dieci & Estep (1991) subsequently gave a broader definition as one in which an ODE 
solver does not detect that the underlying solution is physically unstable. They observed that 
super-stability can occur also when the ODE solver is not super-stable in terms of Dahlquist 
et al. They concluded that the key factor which determines the occurrence of super-stability 
is the iterative solution process for the nonlinear algebraic equations. They indicated that the 
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iterative solution process has its own dynamics, which might be in conflict (as for Newton s 
method) with the dynamics of the problem. They also indicated that super-stability can arise 
because of this fact. Their viewpoint is that the numerical scheme and the methods for solving 
the nonlinear algebraic equations should be considered as a unit. Neglecting this latter aspect, 
and basing step size selection purely on accuracy considerations, leads to faulty analysis. They 
believe that error control strategies for stiff initial value problems ought to be redesigned to take 
into account stability information of the continuous problem. 

In Yee & Sweby (1994, 1995a,b) we exploited the global asymptotic behavior of some of 
the superstable implicit LMMs for constant step sizes. We concentrated on four of the typical 
unconditionally stable implicit LMMs. These are backward Euler, trapezoidal rule, midpoint 
implicit, and three-level backward differentiation (BDF), each with noniterative (linearized 
implicit), simple, Newton and modified Newton iterative procedures for solving the resulting 
nonlinear algebraic equations. A semi-implicit predictor method also was investigated. 

We believe that some of the phenomena observed in our study were not observed by Dieci 
& Estep (1991) or by Iserles (1988). Based on our study, we now give a loose definition of an 
implicit time discretization as having a super-stability property if, within the linearized stability 
limit for a combination of initial data and time step (fixed) or starting time step (using standard 
variable time step control based on accuracy requirements), the scheme stabilizes unstable 
steady states of the governing equations in addition to having the property of Dieci and Estep. 
The definition includes the procedures in solving the nonlinear algebraic equations. This loose 
definition fits the behavior that was observed in Yee & Sweby (1994, 1995a,b), while at the 
same time it fits the framework of dynamics of numerics in time-marching approaches in CFD. 
This is not a re-definition of Dieci and Estep, but rather a clarification of their definition when 
asymptotic numerical solutions of the governing equations are desired. In this case, superstable 
schemes might have the property of a numerical basin of attraction of the true steady state being 
larger than the underlying exact basin of attraction. As can be seen in Yee & Sweby (1994, 
1995a,b), the trapezoidal method is more likely to exhibit this property than the other three 
LMMs. The stabilization of unstable steady states by LMMs was also observed by Salas et al. 
(1986), Embid et al. (1984), Burton & Sweby (1995) and Poliashenko (1995). Section 4.2.2 
gives a summary of the work of Burton and Sweby. 

3.5.2. Implicit LMMs 


The four LMMs and a semi-implicit predictor method considered are 

Implicit Euler Method 


u n+1 = u n + A ts n+1 , 


(3.24) 


u n+1 = u n + Jr(S n + 5 n+1 ), 


Trapezoidal Method 


(3.25) 
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Midpoint Implicit Method 

u n+1 = u" + rS , (3-26) 

3-Level Backward Differentiation Formula (BDF) 

u n+1 = u” + lrS n+1 + ^(« n - (3-27) 

3 o 

Semi-Implicit Predictor Method. 

The semi-implicit predictor method is the same as (3.24) but with an added predictor step using 
the explicit Euler before the implicit step (3.24) to make the final scheme second-order accurate. 
The four methods of solving the resulting nonlinear algebraic equations are as follows. 

Linearization (a noniterative procedure) is achieved by expanding S n+1 as 5 n + 
J(u n )(u" +1 - u n ), where J(u n ) is the Jacobian of 5. 

Simple Iteration is the process where, given a scheme of the form u n+1 = G(u n ,u n+1 ), 
we perform the iteration 

<+.) = < 3 - 28 > 

where u„ +1 = tt n and “(»/)” indicates the iteration index. The iteration is continued either 
until some tolerance between iterates is achieved or a limiting number of iterations has been 
performed. In all of our computations the tolerance “tol” is set as Hu"^ 1 — U (V-\)I I - t0 ^ ^ 
the maximum number of iterations is 15. The major drawback with simple iteration is that for 
guaranteed convergence the iteration must be a contraction; i.e. 

||G(u n ,t>) - G(u n ,to)|| < a||t> — u>||, (3.29) 

where a < 1. Whether or not the iteration is a contraction at the fixed points will influence 
the stability of that fixed point, overriding the stability of the implicit scheme. Away from the 
fixed points the influence will be on the basins of attraction. In other words, “implicit method 
+ simple iteration” behaves like explicit methods. As can be seen in Yee & Sweby (1994, 
1995a,b), numerical results illustrate this limitation in terms of basins of attraction as well. One 
advantage of the “implicit method + simple iteration” over non-LMM explicit methods is that 
spurious steady states cannot occur. Due to this fact, results using simple iteration will not be 
presented here. 

Newton Iteration for the implicit schemes is of the form 


u 


n+l 

(•'+») 


= U 


n+l 

(O 


- F'(u T 


,u 


n+l 1 

(O > 


F(u 


n n+l \ 

»V) h 


(3.30) 


where u£ +1 = u". The differentiation is with respect to the second argument and the scheme 
has been written in the form (for two-level schemes) F(u n ,u n+l ) = 0. 
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Modified Newton iteration is the same as (3.30) except it uses a frozen Jacobian F'(u n , u n ). 
The same tolerance and maximum number of iterations used for the simple iteration are also 
used for the Newton and modified Newton iterations. 

In all of the computations, the starting scheme for the 3-level BDF is the linearized implicit 
Euler. 

We also considered two variable time step control methods. The first one is “implicit Euler 
+ Newton iteration with local truncation error control’’ 


«5J' = 


ti n+1 — ti n+1 4- 


I - 




V = 1 , ... 


(3.31a) 


with 

A t n = 0.9 A< n_I {toll / 1 |u n - - Af"- 1 5(u n )||} 1/2 , (3.31b) 

where the (n + l)th step is rejected if ||u n — u n_1 — A< n-1 5(u n )|| > 2toli. In this case, we 
set A< n = Af n_1 . The value “toll” is a prescribed tolerance and the norm is an infinity norm. 
The second one is the popular “ode23” method 


k 2 = S(u n ) 
k 2 = 5(u" + A< n fci) 
fc, = 5(«" + A< n (Jfe, + jfe 2 )/4) 
u n+1 = u n + A <”(*! +k 2 + 4Jfe s )/6 
Au n+1 = A < n (*! -f k 2 - 2*,)/3 (3.32a) 

with 

A,n=0 - 9Ar ' , /iiS’ (3 - 32b) 

where the (n + l)th step is rejected if ||Au n+1 || > toll max{l, ||« n+1 ||}. In that case, we set 
A t n ~ l = A t n . Again, “toll” is a prescribed tolerance and the norms are infinity norms. We 
also use the “straight Newton” method to obtain the numerical solutions of S(u) = 0, which 
is the one-step Newton iteration of the implicit Euler method of (3.30). 

We mapped out the bifurcation diagrams and numerical basins of attraction for these five 
schemes as a function of the time step for different nonlinear model equations with known 
analytical solutions (scalar and 2x2 systems of autonomous nonlinear ODEs). The computations 
were performed on the CM-5. In general, we preiterated the discretized equations 3000 - 5000 
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steps before plotting the next 3000 iterations. The next section shows some representative 
global asymptotic behavior of these implicit LMMs. 

3.5.3. Numerical Examples 

Selected results in the form of bifurcation diagrams and basins of attraction are shown in 
Figs. 3.10-3.16 for the logistic ODE model (3.7) and Figs. 3.17-3.19 for the second model 
(3.23). In Figs. 3.10, 3.12, 3.14, 3.16-3.19, the left diagrams show the bifurcation diagrams and 
the right diagrams show the bifurcation diagrams with basins of attraction superimposed on the 
same plot. However, Figs. 3.1 1, 3.13 and 3.15 show only the bifurcation diagrams with basins 
of attraction superimposed on the same plot. In all of Figs. 3.10-3.19, the abscissa is r = a At. 

Note that the preselected regions of. Aland the initial data were determined by examining 
a wide range of At and initial data. In most cases, we examined At from close to zero up to 
one million. What is shown in these figures represents some of the Af and initial data ranges 
that are most interesting. Due to this fact, the At and initial data ranges shown are different 
from one method to another for the same model problem. The streaks on some of the plots are 
either due to the non-settling of the solutions within the prescribed number of iterations or the 
existence of small isolated regions of spurious asymptotes. Due to the high cost of computation, 
no further attempts were made to refine their detailed behavior since our purpose was to show 
how, in general, the different numerical methods behave in the context of nonlinear dynamics. 

Due to the method of tracking the basins of attraction, the color coding of the basin of 
attraction associated with a particular asymptote might vary from one vertical line to the next 
vertical line (i.e, from one Af to the next). This is the case for Figs. 3.10, 3.11, 3.14, 3.18 and 
3.19. For example in Fig. 3.10 the basin of attraction (as a function of Af) for the steady state 
ti = 1 is the red region before the appearance of the light blue strip (the first strip). After the 
appearance of the blue strip, it is the region above the curve line that separates the green and 
red regions. (When in doubt, use the bifurcation diagram as a guide and identify the r value 
where the sudden birth of stable spurious asymptotes occurred.) Incidently, for this particular 
discrete map, this envelope and the critical value r which undergoes stabilization of the fixed 
point u = 0 can be obtained analytically. Independently, Arriola (1993) derived the analytical 
form of the envelope. 

The midpoint implicit method behaves in a manner similar to the trapezoidal method. In 
fact their linearized forms are identical. From here on, the midpoint implicit method is not 
discussed. 

As mentioned before, for unconditionally stable LMMs, the scheme should not experience 
any steady bifurcation from the stable branch because unconditionally stable LMMs preserve 
the stability of the stable steady states. This rule does not apply to unstable steady states 
using super-stable LMMs. Before stabilization of the unstable steady state, super-stable LMMs 
typically undergo “inverted period doubling bifurcations” or the “inverted logistic-map-type 
bifurcations” (or crisis in terms of Grebogi et al. 1983). See Fig. 3.10 for —1 < u n < 0.2 for 
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an example. For the ODE (3.23), all of the implicit methods experience at least two inverted 
logistic-map-type bifurcations. From Figs. 3.17-3.19, we can observe that all of the studied 
implicit methods can introduce stable spurious chaos since a logistic-map-type of spurious 
bifurcations occur. 

Figures 3.10, 3.11, 3.14, 3.16-3.19 show other situations where the rest of three implicit 
LMMs and the semi-implicit method stabilize the unstable steady states of the ODEs. It appears 
that the onset of stabilization of the unstable steady states arises in two ways. One way is the 
birth of stable spurious asymptotes or stable spurious chaos in the form of an inverted logistic 
map. The second way is the birth of unstable spurious asymptotes (fixed points other than 
period one) leading to the onset of stabilization of unstable steady states. Although the two 
ways of stabilization of the unstable steady states are similar, their corresponding basins of 
attraction are very different. See Figs. 3.10, 3.14, 3.16-3.19. 

The critical value of A t c for the onset of stabilization is not very large. It is problem 
dependent, method dependent and also procedure dependent (the various ways of solving the 
nonlinear algebraic equations). In most cases, the value of the At c is comparable to or smaller 
than the equivalent of the stability limit of standard explicit Runge-Kutta methods. It is not 
uncommon for the underlying basins of attraction to be larger than the exact basin of attraction 
for At < Af c . 

Among the three procedures, the linearized noniterative forms have a higher tendency to 
stabilize unstable steady states. See Figs. 3.10, 3.14, 3.16-3.19. Here the word “procedure” 
excludes the simple iteration method. Among the three LMMs, the trapezoidal method is the 
least likely to stabilize unstable steady states, but the corresponding basins of attraction can be 
very small and more fragmented than for the other two LMMs. Also, the A t e for stabilization is 
bigger than for the other two LMMs. For a particular LMM not all three procedures for solving 
the nonlinear algebraic equation necessarily stabilize the unstable steady states (see Figs. 3.14 
and 3.15). But, if they do, the pattern or the method for the onset of the stabilization does not 
have to be the same (see Figs. 3.10 and 3.1 1), but the value of A t e is the same for all models 
and methods studied. 

For the case of the semi-implicit predictor method, the onset of the stabilization can be 
accompanied by the birth of other spurious asymptotes (other than steady state). See Fig. 3.16 
for r > 2. It is fascinating to see how complicated the basins of attraction are which compose 
the many disjoint and fractal like regions. Similar behavior is also observed for the “implicit 
Euler + modified Newton” but is less pronounced than for the semi-implicit predictor method. 
See Figs. 3.11 and 3.16. 

Compared with the three implicit LMMs, and independent of the method of solving the 
nonlinear algebraic equation, the semi-implicit method exhibits the smallest basin of attraction 
(compared with the exact basin of attraction of the stable steady state) and is more fragmented 
for At < At c . Aside from the efficiency issue, the implication is that a higher order accuracy 
scheme might not be as desirable for the time-marching approach. 
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Since straight Newton is just a one step “implicit Euler + Newton”, its basins of attraction 
(for At larger than explicit Runge-Kutta) are almost the same, even with more than one step 
iterations. Studies indicated that contrary to popular belief, the initial data using the straight 
Newton method may not have to be close to the exact solution for convergence. Straight Newton 
exhibits stable and unstable spurious asymptotes. Initial data can be reasonably removed from 
the asymptotic values and still be in the basin of attraction. However, the basins can be 
fragmented even though the corresponding exact basins of attraction are single closed domains. 
See Fig. 6.25 of Yee & Sweby (1995a). The cause of nonconvergence may just as readily be due 
to the fact that the numerical basins of attraction are fragmented. If one uses a time step slightly 
bigger than the stability limit of standard explicit methods for the three LMMs, straight Newton 
can have similar or better performance. In fact, using a large At with the linearized implicit 
Euler method or the implicit Euler + Newton procedure has the same chance of obtaining the 
correct steady state as the straight Newton method if the initial data are not known or arbitrary 
initial data are taken. 

Numerical experiments performed on the two variable time step control methods also 
indicated that, although variable time step controls are not foolproof, they might alleviate 
the spurious dynamics most of the time. One shortcoming is that in order to avoid spurious 
dynamics, the required size of At is impractical to use in CFD, especially for the explicit ode23 
method. 

A consequence of all of the observed behavior is that part or all of the flow pattern can 
change topology as the discretized parameter is varied. An implication is that the numerics 
might predict, for example, a nonphysical reattachment flow. Thus even though LMMs preserve 
the same number (but not the same types or stability) of fixed points as the underlying DEs, the 
numerical basins of attraction of LMMs do not coincide with the exact basins of attraction of 
the DEs even for small At. Some of the dynamics of the LMMs observed in our study can be 
used to explain the root of why one cannot achieve the theoretical linearized stability limit of 
the typical implicit LMMs in practice when solving strongly nonlinear DEs (e.g., in CFD). 

Comparing the results with the explicit methods, it was found that aside from exhibiting 
spurious asymptotes, all of the four implicit LMMs can change the type and stability (unstable 
to stable) of the steady states of the differential equations. They also exhibit a drastic distortion 
but less shrinkage of the basin of attraction of the true solution when compared to the standard 
non-LMM explicit methods. Comparing the results of Yee & Sweby (1994, 1995b) with Yee & 
Sweby (1995a), the implication is that unconditionally stable implicit methods are, in general, 
safer to use and have larger numerical basins of attraction than explicit methods. However, one 
cannot use too large a time step since the numerical basins of attraction can be so small and/or 
fragmented that the initial data has to be very close to the exact solution for convergence. This 
knowledge improves the understanding of the basic ingredients needed for a time-marching 
method using constant step size to have a rapid and guaranteed convergence to the correct 
steady state. 
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3.6. Does Error Control Suppress Spuriosity? 

The previous sections discussed mainly the spurious behavior of long time integrations of 
initial value problems of nonlinear ODE solvers for constant step sizes. The use of adaptive 
step size based on local error control for implicit methods was studied by Dieci & Estep (1991). 
Dieci and Estep concluded that for superstable LMMs with local step size error control and 
depending on the procedure for solving the resulting nonlinear algebraic equations, spurious 
behavior can occur. Our preliminary study on the two variable step size control methods 
discussed in the previous section indicated that one shortcoming is that the size of At needed 
to avoid spurious dynamics is impractical (too small) to use, especially for the ode23 method. 
Theoretical studies on the adaptive explicit Runge-Kutta method for long time integration have 
been gaining more attention recently. Recent work by Stuart (1994, 1995), Humphries (1992), 
Higham and Stuart (1995) and Aves et al. (1995) showed that local error control offers benefits 
for long-term computations with certain problems and methods. Aves et al. addressed the heart 
of the question of whether local error control confers global properties of steady states of the 
IVP of autonomous ODEs using adaptive Runge-Kutta type methods. 

Aves et al.’s work is concerned with long term behavior and global quantities of general 
explicit Runge-Kutta methods with step size control for autonomous ODEs. They believed that 
the limit <„ — > oo is more relevant than the limit of the variable step sizes h n — » 0. They 
studied spurious fixed points that persist for arbitrarily small error tolerances r. This type 
of adaptive Runge-Kutta method usually consists of a primary and secondary Runge-Kutta 
methods of different order. Their main result is positive. When standard local error control is 
used, the chance of encountering spuriosity is extremely small. For general systems of ODEs, 
the constraints imposed by the error control criterion make spuriosity extremely unlikely. 
For scalar problems, however, the mechanism by which the algorithm succeeds is indirect — 
spurious fixed points are not removed, but those that exist are forced by the step size selection 
mechanism to be locally repelling (with the relevant eigenvalues behaving like 0(l/r)). 

To be more precise, adaptive time stepping with Runge-Kutta methods involves a pair of 
Runge-Kutta formulae and a tolerance parameter "r", which is usually small. (See the “ode23” 
method (3.32), for example.) Hence a spurious fixed point of the adaptive procedure requires: 

1) A spurious fixed point common to both methods must exist. This is usually easy to achieve 
since the bifurcation diagrams of individual Runge-Kutta methods have so many branches. 

2) This spurious fixed point must be stable. This is much more difficult to achieve - essentially 
since the bifurcation curves for the two methods must intersect tangentially; otherwise there 
will be an eigenvalue of the Jacobian of 0(l/r). 

Aves et al. showed that the probability of 2) occurring is zero. However, for a given pair of 
formulae one can generally construct functions for which it holds (generally stability will only 
hold for r > lower bound). In any event, the basin of attraction of this spurious fixed point will 
only be O(t). These results were derived for scalar problems. 
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In the worst scenario, problem classes exist where, for arbitrary r, stable spurious fixed 
points persist with significant basins of attraction. They derived a technique for constructing 
ODEs for which an adaptive explicit Runge-Kutta method will behave badly. They showed that 
this can be accomplished using a locally piecewise constant function S(ti). Since the disjoint 
pieces can be connected in any manner, S can be made arbitrarily smooth. Hence, smoothness 
of S alone is not sufficient to guarantee that spurious behavior will be eliminated. These 
examples highlight the worst-case behavior of adaptive explicit Runge-Kutta methods. They 
also mentioned that they can construct similar examples involving systems. However, these 
types of examples are somewhat contrived. 

Griffiths is currently working on the application to hyperbolic PDEs. His preliminary results 
(David Griffiths (1996), private communication) showed that it is by no means clear at the 
moment whether stable spurious solutions may be eliminated. The difference is that, unlike 
physical problems governed by nonlinear ODEs, nonlinear PDEs may have wave-like solutions 
rather than fixed points due to the spatial derivatives. 

3.7. A Reaction- Convection Model 

This section further studies the dynamics of selected finite difference methods in the 
framework of a scalar model reaction-convection PDE (LeVeque & Yee 1990) and investigates 
the possible connection of incorrect propagation speeds of discontinuities with the existence 
of some stable spurious steady-state numerical solutions. The effect of spatial as well as time 
discretizations on the existence and stability of spurious steady-state solutions will be discussed. 
This is a summary of the work of Lafon & Yee (1991, 1992). 

A nonlinear reaction-convection model equation in which the exact solution of the governing 
equations are known (LeVeque & Yee 1990) is considered. The model considered in LeVeque 
and Yee is 


tit + au m — oS(u ) 0 < * < 1 (3.33a) 

ti(x, 0) = ti°(*) (3.33b) 

where a and a are parameters, and 5(ti) = — u(l — u)(2 — «). The boundary condition given 
by 


u(0,f) = tio t > 0 (3.33c) 

or, periodic boundary condition given by 

u(0,<) = u(l,t) <>0 (3.33d) 


is used to complete the system. 
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The exact steady-state solutions u* of the continuum PDE (3.33) can be obtained by 
integration by parts of adu* jdx — a5(ti*) which yields 


ax 

— +c = 


a 



= log 


!«*(») - il 1 

W \ u *(.}(2 - **(-))! J ’ 


(3.34) 


where c is the integration constant. 


In the case where the boundary condition is v(0, t ) = u«, there is a unique steady state and 
its value is determined by u 0 . If u 0 is a root of S, (i.e., u 0 = 0, 1, or 2), then the exact steady 
state is constant in * and is equal to u 0 . But if u 0 is not a root of S, then the exact steady state 
satisfies 


x = 



'«*(*) - 1 

/ u 0 (u 0 - 2) 

. («0 - 1) ] 

1 u*(x)(u*(x) - 2). 


(3.35) 


Although, the domain is confined to 0 < at < 1, the steady-state solution is defined for 
0 < * < oo. The limit of «*(*) is 0 as at — ► oo for — oo < u 0 < 0 or 0 < u 0 < 1. The limit of 
«*(*) is 2 if 1 < u 0 < 2, or 2 < u 0 < oo. One can show that this exact steady-state solution is 
stable. 


In the case of the periodic boundary condition where u(0,t) = u(l,f) and u*(0) is not a root 
of S, it can be shown that there exist three exact steady-state solutions; namely «*(*) = 0, 1, 
or 2. One can also show that u* = 0 and u* = 2 are stable while u* = 1 is unstable. 

Denote the basin of attraction for the steady state u* by BA(u*). Then it is obvious that 
IL4(0), the basin of attraction for the steady-state solution u* = 0, is the set of initial data 
«°(*) < 1 for all real values of x. In mathematical notation 

5A(0) = {u° : ti°(at) < 1 V*}. (3.36) 

Similarly, the basin of attraction for the steady state u* = 2 is 

BA( 2) = {u° : «*(*) > 1 V*}. (3.37) 

Later we contrast these exact basins of attraction with the numerical basins of attraction BA( 0) 
and BA( 2) for the various schemes under discussion. 

For the numerical methods, semi-discrete (method of lines) finite difference methods (FDMs) 
and implicit treatment of the source terms (semi-implicit) with noniterative linearization using 
a characteristic form are considered. 


Spatial^ Discretization tj Let uj(t) represent an approximation to «(*,-, t ) where xj = j Ax and 
j = 0, ..., J with Ax = 1/J the uniform grid spacing ( J + 1 grid points). Let the parameter 

a 

a Ax 


Pi = 


(3.38) 
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Define the flow residence time in a cell r c = Ax/a (characteristic time due to convection) 
and the time required by the reaction to reach equilibrium r r = 1/a. Then a simple physical 
interpretation of the parameter pi comes from the fact that pi is equal to the ratio t v /t c . Note 
also that this ratio is the inverse of a Damkohler number. Therefore, the parameter pi is a 
measure of the stiffness of the problem. The smaller pi is, the stiffer the problem becomes. 

A semi-discrete approximation (for a chosen spatial discretization for the convection and 
source term) of the reaction-convection PDE (3.33a) is then 

-—=F, (3.39) 

a dt 

where U is the vector whose components are u,(f), 1 < j < J- The function F is a 
discrete J-dimensional vector which .depends on the grid function U, the parameter p x , and 
the particular spatial finite difference approximation. For simplicity the commonly used spatial 
discretizations (the first-order upwind (UP1), second-order upwind (UP2) and second-order 
central (C2) schemes) are considered for the convection term, and the pointwise evaluation 
(PW), upwind interpolation (UI), and mean average between two neighboring grid points (MA) 
are considered for the (spatial) numerical treatments of the source term. The combination of 
the three numerical treatments of the source term and the three basic schemes used for the 
discretization of the convection term yields 9 spatial finite difference approximations for the 
reaction-convection PDE (3.33a). 

The expressions of the elements fj of F corresponding to the 9 spatial difference approxima- 
tions for (3.33a) and (3.33c) are given below, where we use the obvious notations u_i = uj_ x , 
u 0 = uj and ui = u J+ i for the periodic boundary condition (3.33c). 

1. First-order upwind for convection - pointwise evaluation for source term (UP1PW) 

fj(U) = -Pi(uj - uj-i) + S(uj). (3.40a) 

2. Second-order upwind for convection - pointwise evaluation for source term (UP2PW) 

fi(U) = -Pi(^tii - 1 + \vj-t) + £(«;). (3.40b) 

3. Second-order central for convection - pointwise evaluation for source term (C2PW) 

fj(l o = — 2 Pi( u i+i “ u i-0 + £(“>)• (3.40c) 

4. First-order upwind for convection - upwind evaluation for source term (UP1UI) 

fj(U) = — pi(«j — u i-i) + 0S(uj-i) + (1 — 0)5(tij). (3.41a) 

5. Second-order upwind for convection - upwind evaluation for source term (UP2UI) 

fj(U) = — Pi(^«j - u j-i + 4 U i-*) + ^( u i-i) + (1 - 0)S( u j)‘ 


(3.41b) 
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6. Second-order central for convection - upwind evaluation for source term (C2UI) 

(3.41c) 

7. First-order upwind for convection - mean average evaluation for source term (UP IMA) 

fj( u ) = -Pi(«i “ u i - 0 + \ s ( u j+ 1) + s ( u i- 0 • (3.42a) 

8. Second-order upwind for convection - mean average evaluation for source term (UP2MA) 

fj( u ) - ~ 1 + + \ 5 ( w i-0 + 5 ( u i+i) • (3.42b) 

9. Second-order central for convection - mean average evaluation for source term (C2MA) 

fj ( u ) = ~\pi{ u i+i ~ u j-i) + \ s ( u j- 1) + ^(«i+i) • (3.42c) 

The parameter 9 associated with the upwind interpolation of the source term in formula (3.41) 
is an extra parameter in the discretization, lying between 0 and 1. For ease of referencing, 
the above 9 spatial discretizations will be denoted by the symbols UP1PW, UP1UI, UP1MA, 
UP2PW, UP2UI, UP2MA, C2PW, C2UI and C2MA. 

Time Discretizationt : To contrast the nonlinear behavior between LMMs and explicit Runge- 
Kutta methods, for simplicity, the explicit Euler and the modified Euler schemes are considered. 

Let u” represent an approximation of u(j A®, nAt) with a fixed time step Af. Also let U n 
denote a vector with elements u”. Then the fully discrete approximation of the PDE (3.33a) 
with explicit Euler time discretization is 

V n+i = u n +p 2 F(U n ), ( 3 . 43 ) 

where the parameter pj is simply related to the time step through 

P 2 = aAt. (3-44) 

With modified Euler time discretization, the fully discrete approximation is 

U n+1 = U n +P2F(U) ; U = U n + ^ F(U n ). (3.45) 

In the following the above 18 fully discrete approximations will be denoted for ease of 
reference by the symbols UP1PW/EE, UP1PW/ME, etc. (where EE stands for explicit Euler 
and ME for modified Euler). In all of the computations 0 = c = in (3.41) is used. 
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FDMt Bated on the Characteristic Form : A more physical approach to updating the grid 
value uj at time level n + 1 is to trace back the characteristic passing through (xj, (n + l)Af). 
Denote by (i,nAi) the coordinates of the point on the characteristic at time nAt. Along this 
characteristic, the problem reduces to solving the ODE 

v T = aS'(v) nAt < r < (n + 1)A t (3.46a) 

with the initial condition 


t>(nAf) = u n (£), (3.46b) 

where u"(i) denotes the value of u at point x and time level nAt, obtained by some 
interpolation method on the adjacent grid function u". For this approach, explicit and implicit 
time discretizations of (3.46a) or, equivalently, explicit and implicit treatments of the source 
term of the PDE (3.33a) are considered. To facilitate the analysis, the convection part of the 
PDE (3.33a) is handled in an explicit way. This means that for the homogeneous part of the 
PDE (3.33a) all the underlying schemes are under the CFL restriction 

c=~=p in < 1. (3.47) 

Ax 

An immediate consequence is that x lies between Xj-i = (j — 1)A* and xj = j Ax. Then, 
from a linear interpolation we get 

u n = u n (x) = cu^_ 1 +(l-c)u^. (3.48) 

For fully explicit schemes, the same two explicit time discretizations for the method of lines 
approach are considered here. With explicit Euler, the fully discrete approximation (denoted by 
CHA/EE) takes the form 


u? +1 = cu7_, + (1 - c)«7 +p 2 5(« n ). (3.49) 

With the modified Euler time discretization the fully discrete approximation (denoted by 
CHA/ME) takes the form 


u? +1 = c«7_ 1 + (1 - c)u? + piSi*), (3.50) 

where 

0 = c«7. 1+ (l-c)u; + ^S(fi"). (3.51) 

For a less restricted time step, the implicit Euler (IE) and the trapezoidal implicit scheme 
(T) are considered for the source term. With implicit Euler, the fully discrete approximation, 
denoted by CHA/1E, takes the form 

u] + ' -p J S(u7 +1 ) = u", 


(3.52) 
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while with the trapezoidal rule, the fully discrete approximation, denoted by CHA/T, takes the 
form 


u” + ' - f 5(u" +1 ) = + (1 - + T s ( 4 ”)- 


(3.53) 

2 - \ ~j / j— a \ '^ 2 V/ 

The corresponding linearized implicit scheme CHA/LIE associated with the fully discrete 
approximation CHA/IE is given by 


u 


r»+l 


= «7 + 


cuj_i -cu? +P25 (u?) 
l-ptS'iu'J) 


(3.54) 


with S’ — dSjdu — -2 + 6u - 3u 2 . The linearized trapezoidal method CHA/LT is 

cu^_ x - cu? + ^[5(i»?) + 5(fi w )] 


+ 




(3.55) 


3.7.1. Spurious Asymptotes of Full Discretizations 

Besides the three exact steady-state solutions, depending on the numerical methods, either 
the spatial discretizations and/or the time differencing can independently introduce spurious 
asymptotic numerical solutions (see Lafon & Yee (1991) for a detailed analysis). Bifurcation 
diagrams and stability analysis for the exact and spurious asymptotes of the above schemes and 
source term treatments were discussed in Lafon & Yee (1991, 1992). Interested readers are 
referred to these references for more details. 

Since the explicit Euler and the implicit methods are LMMs, no spurious steady state due 
to time discretizations are possible. But, consider for example the various FDMs involving 
modified Euler time differencing (method of lines or characteristic form) and look for the simple 
case of spatially invariant steady states (i.e., the value of Uj is the same for all j, 1 < j < J). 
Then it is found that for such FDMs, the value u* of a spatially invariant steady-state must 
satisfy 

«* + ^5(u*) = u, 5(u) = 0. (3.56) 

A 

It can be shown that (3.56) admits the following 9 solutions 


3± Jl + — 

, 1, 1±4 

1 2 1 

,/l H 1 2, - 

!±i 

h- 

V Pi J 


If Pi 2 


If Pi J 


in which u * = 0, 1 or 2 are the exact steady-state solutions while the rest are spurious 
steady-state numerical solutions introduced by the modified Euler time discretization. 
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3.7.2. Linearized Behavior vs. Nonlinear Behavior 

To illustrate the differences between the linearized analysis and the nonlinear solution 
behavior. Fig. 3.20 shows the spectral radius around the exact steady-state solution u* = 2 
and the bifurcation diagram obtained with initial data Z7 # = (2. ,2.1, 2., 2. 1,2. ,2.1, 2.) for the 
scheme UP1UI/EE and pi = 7. Similar results for the scheme UP1UI/ME are shown in 
Fig. 3.21. For pi = 7, the scheme UP1UI/ME exhibited two disjoint linearized stability 
intervals. From Fig. 3.21 we observe that outside these stability intervals the scheme 
does not necessarily diverge (as indicated by the linearized analysis) but can converge to 
spurious asymptotic solutions. The spurious behavior is very different between the two time 
discretizations employing the same spatial discretization and initial data. The modified Euler 
exhibited spurious steady states, whereas the explicit Euler exhibited spurious asymptotes other 
than spurious steady states. 

Another example which shows the importance of nonlinear analysis is that of two schemes 
that exhibit the same linearized behavior yet have different nonlinear behavior (true behavior). 
For example, the linearized stability analyses for schemes UP1UI/EE and CHA/EE (see Figs. 
5a and 7a of Lafon & Yee (1992)) are identical even though the bifurcation diagrams obtained 
with the same initial data and pi = 7 are different (see Figs 19a and 20a of Lafon & Yee 
(1992)). 


3.7.3. Spurious Steady States and Nonphysical Wave Speeds 

The possible connection of the numerical phenomenon of incorrect propagation speeds of 
discontinuities with the existence of some stable spurious steady states introduced by the spatial 
discretization is discussed here. The boundary condition (3.33c) and the following piecewise 
constant initial data 


u 


o 



2 

0 


0 < x < Xd, 

Xd < x < 1 


(3.58) 


are considered. The constant value denotes the location of the discontinuity. The exact 
solution of (3.33a,b,c) with initial data (3.58) is simply 


u(as, t) = ti°(* — at) (3.59) 

which is a wave traveling at constant speed a. 

For an explanation of how numerical methods applied to this model PDE are likely to produce 
wrong speeds of propagation for the initial data (3.58), the reader is referred to LeVeque & 
Yee (1990). This phenomenon is due to the smearing of the discontinuity caused by the spatial 
discretization of the advection term. This introduces a nonequilibrium state into the calculation. 
Thus, as soon as a nonequilibrium value is introduced in this manner, the source term turns on 
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and immediately tries to restore equilibrium, while at the same time shifting the discontinuity 
to a cell boundary. 

For simplicity, consider the first-order upwind spatial discretization with the explicit Euler 
time discretization for (3.33a,c) and (3.58). Assume an equal spatial increment of J intervals so 
that the discretized initial data associated with (3.59) is 


| 2 1 < j < K 


(3.60) 


with the index K depending on the constant *<j. In addition, assume that the convection speed 
a — 1 so that At = c/J. In the computation J — 20 and the Courant number c = .05. With 
these assumptions, only the stiffness a of the source term is a free parameter. Define the average 
wave speed W for the numerical solution as follows 


W = 


A* 

2nAf 


pJ+i 


L i= i 


J + 1 


Vi 


j = 1 


(3.61) 


where the average is taken on the time interval 0 to nAt. Figure 3.22 shows this average speed 
versus pi (proportional to 1/a). It reveals the fact that when pi becomes large (or, similarly 
when the source term is not stiff) the numerical wave speed tends to the exact solution (a = 1), 
while for pi < .25 the numerical solution is a standing wave (the average speed being 0). 


This zero wave speed is indeed a stable spurious steady-state solution of the discretized 
equation, but not a solution of the continuum PDE. Since the explicit Euler time differencing 
has the property of not producing spurious steady-state numerical solutions, this spurious stable 
steady-state numerical solution must have been introduced by the spatial discretization. In 
fact, it is evident from the bifurcation diagram shown in Fig. 2 of Lafon & Yee (1991) 
that this spurious steady state lies on the stable branch originating at pi =0 from the state 

til — 2 } ....j vjf = 2, ujif + i — 0,....,t ij — 0. 


3.7.4. Numerical Basins of Attraction 


In order to show the global nonlinear solution behavior of the schemes, numerical basins 
of attraction (of the underlying schemes) are compared with each other as well as with the 
exact basin of attraction for u* = 2 (BA( 2)). Due to the complexity and CPU intensive nature 
of the computation, unlike the detailed basins of attraction presented in Yee & Sweby (1994, 
1995a,b), only a fraction of the basins of attraction are computed. 

To compute a partial view of some numerical basins of attraction for u* = 2 (BA( 2)), a set 
of initial data in a two-dimensional plane was selected. Even with this restriction, the analysis 
is still very complex and the computation is CPU intensive; it was performed only for the case 
J = 4 (5 grid points). The equation of the chosen plane is 


Uj — 2 j u 4 — 2. 


(3.62) 
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A set of 121 initial data in the plane (3.62) were obtained by confining u 2 and u s in the square 
1.1 < u 2 < 3.1, 1.1 < tij < 3.1 with a uniform increment of Au 2 = Auj = 0.2. For each 
initial datum, 5000 preiterations were performed. The asymptotic behavior is plotted in the 
plane. Results are shown in Figs. 3.23 - 3.25 In all of the figures, open triangles 
belong to the numerical basin of attraction BA( 2). Dark squares belong to the numerical basin 
of attraction of a spurious steady state. Dark circles are the numerical basins of attraction of 
other (exact or spurious) asymptotic solutions, while a blank space denotes a divergent solution. 
Note that the whole square domain 1.1 < u 2 < 3.1, 1.1 < u 2 < 3.1 is inside the true basin of 
attraction of the exact steady state u* = 2 and therefore the true behavior should be an open 
triangle for all the initial data considered. The study in Lafon & Yee (1991) indicated that 
the modified Euler time differencing has a smaller attractive region than the explicit Euler for 
the domain considered even though in terms of linearized stability and accuracy, the modified 
Euler exhibits an advantage over the explicit Euler. 

For example, for pj = 0.1 and p 2 = 0.5, implicit treatments of the source term exhibit a 
larger attractive region of initial data than the explicit treatments of the source term. However, 
as the time step is further increased, an adverse behavior is observed contrary to the common 
belief that implicit schemes can operate with much higher time step and still produce the 
correct steady-state numerical solution. Since the source term is handled implicitly, the classical 
guideline for the time step constraint is given by the explicit discretization of the convective term 
(homogeneous part of the PDE (3.33a)), or equivalently by the CFL constraint c = P 1 P 2 < 1. 
Thus, the permissible time step for the implicit treatments of the source term is larger than 
explicit treatments of the source term. However, it is evident from the computation shown on 
Figs. 3.23 - 3.25 for implicit schemes CHA/LIE and CHA/LT with a Courant number set equal 
to 0.3 (p 2 = 3 ) and 1 (j >2 = 10 ) respectively, that these implicit schemes no longer give the 
correct asymptotic solution, in particular for the scheme CHA/LT. 

3.8. Time-Accurate Computations 

In the examples chosen by Lorenz (1989), he showed that numerical chaos always precedes 
divergence of a computational scheme. He suggested that computational chaos is a prelude 
to computational instability. Poliashenko & Aidun (1995) showed that this is not a universal 
scenario. In previous sections, we have shown that, depending on the initial data, time steps and 
grid spacings, numerics can introduce spurious asymptotes and chaos. Using a simple example, 
Corless (1994b) showed that numerics can suppress chaotic solutions. The work of Poliashenko 
and Aidun also discussed spurious numerics in transient computations. Adams et al. (1993) 
discussed spurious chaotic phenomena in astrophysics and celestial mechanics. They showed 
that the source of certain observed chaotic numerical solutions might be attributed to round-off 
errors. Adams also discussed the use of interval arithmetics (interval mathematics or enclosure 
methods) to avoid this type of spurious behavior. Moore et al. (1990) discussed the reliability 
of numerical experiments in thermosolutal convection. Keener (1987) discussed the uses and 
abuses of numerical methods in cardiology. 
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It is a common misconception that inaccuracy in long-time behavior of numerical schemes 
poses no consequences for transient or time-accurate solutions. This is not the case when one is 
dealing with genuinely nonlinear DEs (Jackson 1989). For genuinely nonlinear problems, due 
to the possible existence of spurious solutions, larger numerical errors can be introduced by the 
numerical methods than one can expect from a local linearized analysis or weakly nonlinear 
behavior. Lafon & Yee (1991) illustrated this connection. Section 3.7.3 summarizes their 
result The implication is that it might be possible that a continuum consists of no steady state 
but the numerics might induce spurious equilibrium states or other asymptotes. With certain 
combinations of initial data, time step and grid spacing, the time-accurate computation can 
get trapped in a spurious standing wave or spurious time-periodic solution. Sections 5.4-5.6 
illustrate three different types of spurious behavior in unsteady CFD computations (Yee et al. 
1997). 

The situation can become more intensified if the initial data of the DE is in the basin of 
attraction of a chaotic transient (Grebogi et al. 1983) of the discretized counterpart. In fact, it 
is possible that part or all of the solution trajectory is erroneous. Section V shows a practical 
example of a chaotic transient near the onset of turbulence in direct numerical simulations of 
channel flow by Keefe (1996). Since numerics can introduce and suppress chaos, and can also 
introduce chaotic transients, it casts some doubts in relying on numerical tests for the onset of 
turbulence and chaos. 

There has been much debate on the overall accuracy away from shocks of high-resolution 
shock-capturing methods. Donat (1994) addressed this issue from a theoretical standpoint while 
Casper & Carpenter (1995) illustrated it with a shock induced sound wave model. Casper and 
Carpenter concluded that there is only first-order accuracy downstream of the sound-shock 
interaction using a spatially 4th-order ENO scheme. Sections 5.2 and 5.3 illustrate two additional 
types of spurious numerics in transient computations. 


IV. Spurious Dynamics in Steady-State Computations 

Any CFD practioner would agree that making a time-marching CFD computer code converge 
efficiently to a correct steady state for poorly understood new physical problems is still an art 
rather than a science. Usually, after tuning the code, one still encounters problems such as 
blow-up, nonconvergence, nonphysical nature, or slow convergence of the numerical solution. 
Some of these phenomena have been reported in conference proceedings and reference journals, 
but the majority have been left unreported. Although these behaviors might be caused by factors 
such as poor grid quality, under-resolved grids, improper numerical boundary conditions, etc., 
most often they can be overcome by employing standard procedures such as using physical 
guidelines, grid refinement, improved numerical boundary treatments, halving of the time step, 
and using more than one scheme to double check if the numerical solution is accurate and 
physically correct. However, these standard practices alone may sometimes be misleading, not 
possible (e.g., too CPU intensive) or inconclusive due to the various numerical uncertainties 
(see section I) that can be attributed to the overall solution process. Consequently, isolation of 
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the sources of numerical uncertainties is of fundamental importance. Section III isolates some 
of the spurious numerics for elementary models. Complementing the phenomena observed in 
Section III, this section illustrates examples from CFD computations (Yee & Sweby 1996a). 
We concentrate mainly on the convergence issues that are contributed to the spurious dynamics 
that are inherent in the schemes. Section 4.2. 1 was written by Bjorn Sjogreen of Royal Institute 
of Technology, Sweden. 

4.1. A 1-D Chemically Relaxed Nonequilibrium Flow Model 4 

This section discusses the analysis of numerical basins of attraction for the simulation of a 
1-D chemically relaxed nonequilibrium flow model for a [N 2i N) mixture. Sweby et al. (1995), 
Yee & Sweby (1996a) and Yee et al. (1997) studied the spurious behavior of this model for 
six different explicit Runge-Kutta methods and a semi-implicit form of the implicit Euler and 
trapezoidal methods. This type of flow is encountered in various physical situations such as 
shock tube experiments (the mixture behind the shock being in a highly nonequilibrium state) or 
a high enthalpy hypersonic wind tunnel. Under these assumptions the model can be expressed 
as a single ODE, 

^ = 5(p,T,z), (4-1) 

where z is the mass fraction of the N 2 species, p is the density of the mixture, T is the 
temperature and * is the 1-D spatial variable. There are two algebraic equations for p and 
T. This system consists of a large disparity in the range of parameter values (many orders of 
magnitude) and is stiff and highly nonlinear. 

The derivation of the model is as follows. The one-dimensional steady Euler equations for a 
reacting (JV 2 , N) mixture are 



(4.2a) 

e 

II 

(4.2b) 

o 

II 

+ 

(4.2c) 

— u(E + p) — 0, 
dx 

(4.2d) 

where (4.2a) is the balance equation for the N 2 species and w Ni is the production rate of the 
N 2 species with density p Nl . The variables p, u, E and p are density, velocity, total internal 
energy per unit volume, and pressure, respectively. 


4 We would like to acknowledge Andre Lafon for the original formulation and the earlier study used 
in this section; presented at the ICFD Conference on Numerical Methods for Fluid Dynamics, April 
5-6, 1995, Oxford, UK. 
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The production rate wpf a of species N 2 is the sum of the production rates for the two reactions 


JV 2 + N 2 2JV + N 2 (4.3a) 

N 2 + N ^ 2N + N (4.3b) 

and is computed using Park’s model (Park 1985) that has been used extensively for hypersonic 
computations. See the Workshop on Hypersonics (1991) for some discussion. These reaction 
rates involve an equilibrium constant, K eq (see below), which is determined by a polynomial 
fitting to experimental data, and as such is only valid for a certain range of temperatures. In 
particular, a cut-off value has to be introduced for low temperatures, a typical choice being 
T min = 1000 K (Mulard & Moules 1991). 

The systems (4.2) and (4.3) must be closed by a thermodynamic representation of the 
mixture. Here a simple model with no vibrational effects has been chosen. The details have 
been omitted for brevity. 

Equations (4.2b)-(4.2d) simply integrate to give 

= goo, (4.4a) 

pu 2 +p = P oot (4.4b) 

H = 111 = Hoot ( 4 . 4c ) 

P 

where H is the total enthalpy and goo, -Poo and JT*, are all constants. Finally, denoting the mass 
fraction of the N 2 species by 


z = 


PN» 

P 


(4.5) 


and using Park’s reaction rate model and the thermodynamic closure, (4.1) can be written as 


dz 

dx 


S(p,T,z) 
1 

M 1 q ao 


2 rpB 

p 1 exp 


(- 1 ) 


X 


\a.A\z(l — zf — A 2 z 2 + 2aA 2 (l — z) s — 2A 2 z(l — z)], 


(4.6a) 


where 

a = — Cr — , = 10® exp [cj + c 2 Z + cjZ a + c^Z* + c s Z 4 ], Z = —— . (4.6b) 
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The density p is obtained from 


> (8-2*)0j-(lO-Sz)P.0) +2(2-r)[ff_-(l-z)«5] =» (4.6c) 


and the temperature T from 


and the pressure from 


The model uses the constants 


T = 


M ip 

R( 2 — z)p 



(4.6d) 

(4.6e) 


■ a = 3.898 1 rAi = 3.7 x 10 18 ' 

Cj = -12.611 A 2 = 1.11 x 10“ 

c, = 0.683 B = -1.6 

c 4 = -0.118 9 - 1.132 x 10 8 

c 6 = 0.006 e\ = 3.355 x 10 7 

.Mi = 28x10* J IR = 8.3143 . 

The input parameters goo, Poo and IZoo are set equal to 0.0561, 158,000 and 27,400,000, 
respectively. A limitation of the model is T > T m ; w = 1000JT. The acceptable root of (4.6c) is 
taken to be real and positive. In addition, solutions are nonphysical if z g [0, 1], if p < 0 or if p 
is complex. 

In the integration of (4.1), the spatial variable * acts as a time-like variable. The asymptotic 
state is the equilibrium state given by 5(p, T,z) = 0. Equation (4.6a) was integrated using 
the Euler, modified Euler (R-K 2), improved Euler (R-K 2), Heun (R-K 3), Kutta (R-K 3) and 
4th-order Runge-Kutta (R-K 4) schemes, and a semi-implicit version of the implicit Euler and 
trapezoidal methods. 

There are two strategies possible when implementing these schemes. One is to freeze the 
values of p and T at the beginning of each step when calculating S(p,T,z) at the intermediate 
stages. The other is to update the values at each evaluation of the function S. The results 
presented here employ the latter strategy since this is the more proper implementation; however, 
it is interesting to note (see below) that results obtained by freezing p and T for intermediate 
calculations exhibit a slightly richer dynamical structure. Due to the complexity of the equation 
and the coupling of the unknowns, the implicit method was implemented by treating the z 
variable implicitly and the rest of the unknowns explicitly. 

In each case, the computations were performed for a range of initial z and integration steps 
A*. For each fixed Ax and each initial datum, the discretized equations were preiterated 



52 


1000 steps before a full bifurcation diagram (of the asymptotic states) together with basins of 
attraction were produced. The preiterations are necessary in order for the solutions to settle to 
their asymptotic values. To obtain a bifurcation diagram with numerical basins of attraction 
superimposed, the preselected domain of initial data and the preselected range of the A* 
parameter are divided into 256 or 512 equal increments. We keep track where each initial datum 
asymptotically approaches, and color code each basin according to the individual asymptotes. 
Figures 4. 1 and 4.2 show the results obtained from these computations. Due to the fact that 
for each A*, only two distinct basins of attraction are present for all of the computations, only 
the grey-scale version of these plots are shown. In all of these plots the shaded region denotes 
the basin of attraction in which combinations of “initial” upstream input z values and step 
size Ax converge to the stable asymptotes of the discretized equations, depicted by the solid 
black line or black dots. The unshaded regions indicate regions of upstream initial input where 
the combinations of upstream input z and Ax do not converge or converge to a nonphysical 
solution of the problem (see condition below (4.6e)). As can be seen in all cases, there is a 
drastic reduction in the basins of attraction with just a slight increase in the grid spacing. (The 
axis scale is 10 -5 !) Note that the allowable upstream initial input (exact basin of attraction) for 
the governing equation (4. 1 ) is 0 < z < 1. 

Explicit Runge-Kutta Methods : The explicit Euler scheme (Fig. 4.1a) obtains the 
correct equilibrium state up to its linearized stability limit, where there is a very small region of 
period 2 spurious solutions before it diverges. Similar behavior is observed for the improved 
Euler (Fig. 4.1c) and Kutta (Fig. 4.1e) schemes, the latter also exhibiting a much more 
constricted basin of attraction for any given A*. The Heun scheme (Fig. 4. Id) exhibits a distinct 
region where stable spurious periodic solutions occurred just above the linearized stability limit. 

As is typical with the modified Euler scheme (Fig. 4.1b), a transcritical bifurcation occurs 
at its stability limit which leads to a spurious (A* dependent) solution near the stability limit. 
Note also the solid line at about z = 0.25 down on the plot, outside of the shaded region. 
This appears to be an unstable feature picked up by our method of asymptotic equilibrium state 
detection (comparison of initial data with the 1000th iterate) and is unlikely to arise in practical 
calculations unless the initial data are on this curve. The R-K 4 scheme (Fig. 4. If) also exhibits 
a transcritical bifurcation at the linearized stability limit; however, this is discernible more by 
the sudden narrowing of the basin of attraction since the spurious asymptotic state varies only 
slightly with A*. 

If the values of p and T are frozen for intermediate calculations, the dynamics are somewhat 
modified. All schemes with the exception of explicit Euler have a slightly larger basin of 
attraction for values of A* within the stability limit and all schemes have period two behavior at 
the stability limit, there being no transcritical bifurcations for any of the schemes. The modified 
Euler scheme also has embedded period doubling and chaotic behavior below the linearized 
stability limit. 

Semi-implicit Methods : Due to the complexity of the equation and the coupling of the 
unknowns, the implicit method was implemented by treating the z variable implicitly and 
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the rest of the unknowns explicitly. Only the linearized version (noniterative) of the implicit 
Euler and trapezoidal formula is considered. To aid in the discussion, bifurcation diagrams 
and the bifurcation diagrams with the basins of attraction superimposed on the same plot are 
shown in Fig. 4.2 using this semi-implicit approach. Both methods behave in a manner 
similar to their explicit counterparts, except that the reduction in the basins of attraction (as 
Ax increases) is considerably less than the explicit Runge-Kutta methods for the implicit 
Euler but only slightly larger for the trapezoidal method. In addition, the semi-implicit Euler 
approach exhibits a richer dynamical structure than the rest of the methods studied. Besides a 
definite region where period two spurious solutions bifurcate from the true branch of the steady 
state solutions near A* = 3.1 x 10 “ 5 , spurious periodic solutions up to period 480 occur for 
smaller A* (see the two dark regions for 2.34 x 10 -5 < A* < 2.7 x 10" 5 ). These spurious 
asymptotes appear to be chaotic-like. In other words, it is possible that for the same grid 
spacing but different initial input, the numerical solution converges to two different solutions. 
For 2.34 x 10 -5 < Ax < 2.7 x 10 -5 , one of the numerical solutions is spurious (initial data 
reside on the shaded region or on the two darker regions). For A* > 3.1 x 10 -5 , all of the 
numerical solutions are spurious. 

The above computations illustrate the sensitivity of the allowable upstream initial inputs to 
the slight increase in the grid spacing. In other words, with a slight increase in the grid spacing, 
the allowable upstream initial inputs quickly become “numerically unphysical’’. Although the 
dynamical behavior of the studied schemes is perhaps not as rich as in some of simple examples 
discussed in Yee et al. (1991), Yee & Sweby (1994, 1995a,b), and Lafon & Yee (1996a), 
spurious features can still occur in practical calculations and so care must be taken in both 
computation and in interpretation. 

4.2. Convergence Rate & Spurious Dynamics of High-Resolution Shock-Capturing Schemes 

We have seen in Section III elementary examples and references cited therein on how the 
proper choice of initial data and the step size combination can avoid spurious dynamics. Yet 
for other combinations the numerical solutions can get trapped in a spurious limit cycle. We 
have also seen that the convergence rates of the schemes are greatly affected by the step sizes 
that are near bifurcation points. Here we include the dynamics of full discretization of two 
nonlinear PDE examples. The spatial discretizations are of the high-resolution shock-capturing 
type (nonlinear schemes). This includes TVD and ENO schemes. Section 4.2.1 discusses how 
this nonlinear scheme affects the convergence rate of systems of hyperbolic conservation laws. 
Section 4.2.2 illustrates the existence of spurious asymptotes due to the various flux limiters 
that are built into TVD schemes. 


4.2.1. Convergence Rate for Systems of Hyperbolic Conservation Laws 

This section summarizes the results of Engquist & Sjogreen (1995) and Sjogreen (1996, 
private communication). These results concern the convergence rate for discontinuous solutions 
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of a system of nonlinear hyperbolic conservation laws. For a scalar nonlinear conservation law, 
the characteristics point into the shock. According to the linear theory of Kreiss & Lundqvist 
(1968), dissipative schemes damp out errors propagating backwards against the direction of the 
characteristics. Thus it is reasonable to expect that the locally large errors at the shock stay in 
a layer near the shock. In numerical experiments we usually obtain 0(h p ) convergence away 
from the shock with difference schemes of formally p th order. 

For the systems case, the same reasoning from the scalar conservation law cannot be applied. 
In this case, we have other families of characteristics intersecting the shock causing the situation 
to be more involved. Thus it is possible that the large error near the shock propagates out into 
the entire post shock region by following a characteristic which emerges from the shock. 

This effect cannot be seen in a simple scalar Riemann problem (problem with jump initial 
data), because exact global conservation determines the post shock states. The system model 
problem, taken from Engquist & Sjogreen (1995), 


u t + (u 2 /4), = 0, — oo < x < oo, 0 < t (4.7a) 

v t + Vm + 3(u) = o, y(u) = (ti + l)(u - l)(l/2 - u) (4.7b) 

gives an example of propagation of large errors. The function g(u) has the properties 
<j(l) = fl(-l) = ff(l/2) = 0, and j(u) ^ 0 for -1 < u < 1 except at u = 1/2. The initial 
data was given as 


r i x < o 

uo(z) = | ft » ®*(*) = 1 (4.7c) 

( — 1 x > 0 

so that the exact solution of the u equation is a steady shock. The eigenvalues of the Jacobian 
matrix of the flux (u 2 / 4,v) T for (4.7) are Ai = u/2 and A a = 1 . The eigenvalue Ai = u/2 
corresponds to a strictly nonlinear field, and A 2 = 1 corresponds to a linearly degenerate field. 

With this initial data (4.7c), it gives rise to a steady 1-shock, with the 1 -characteristics having 
a slope 1/2 to the left of the shock and a slope -1/2 to the right of the shock. These thus intersect 
the shock when time increases. The 2-characteristics of the linear field have slope 1 on both 
sides of the shock. These characteristics thus enter the shock from the left and exit to the right. 
The t> component of the solution is passively advected along the 2-characteristics. When these 
characteristics exit from the shock at * = 0, an error, coming from poor accuracy locally at the 
shock, is picked up and advected along with the solution into the domain x > 0. The shock 
curve * = 0 ( in the x-t plane ) acts as an inflow boundary for the domain * > 0. The error 
coming from the shock is similar to an error in given inflow data, and is therefore not affected 
by the numerical method used in the interior of the domain. Thus is not surprising that this 
error is of first order, even when the equation is solved by a method of higher formal order of 
accuracy. 
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Figure 4.3 shows the numerical solution, computed by a second-order accurate ENO method 
using 50 grid points at the time T = 5.68. The points in the shock give a large error which is 
coupled to the u equation through g(u). The exact solution for v is 1. Numerical investigation 
of the convergence rate of the error in v to the right gave the exponent 1.047. One thus has 
first-order convergence for this second-order accurate method. 

Similar effects can be seen in computing the quasi one-dimensional nozzle flow. Sjogreen 
(1996) computed the solution on the domain 0 < * < 10 for a nozzle with the following cross 
sectional area variation 


A(x ) = 1.398 + 0.347 tanh(0.8* - 4). (4.8) 

This problem is studied in Yee et al. (1985) for a class of explicit and implicit TVD schemes. 
The solution has a steady shock in the middle of the domain. Figure 4.4 shows the error in 
momentum for the steady-state solution on grids of 50, 100, and 200 points for a fourth-order 
ENO scheme and a second-order TVD scheme. For the fourth-order method, the convergence 
exponent is 3.9 before the shock and 1.0 after the shock, when going from 100 to 200 points. 
For the second-order TVD the same quantities have the values 2.2 and 1.1 respectively. 

Sjogreen (1996) recently conducted the same numerical study for the two-dimensional 
compressible Euler equations for a supersonic flow past a disk with Mach number 3. The 
equations were discretized by a second-order accurate uniformly nonoscillatory (UNO) scheme 
(Harten 1986), which unlike TVD schemes, is formally second-order everywhere including 
smooth extrema. He computes the error in entropy along the stagnation line for the steady-state 
solution on grids with 33 x 17, 65 x33, and 129 x65 grid points. The result is shown in Fig. 4.5, 
where the error and convergence exponent in the region behind the bow shock are plotted. The 
convergence exponent is between the 65 x 33 and 129 x 65 grids. The disk has radius 0.5, and 
it is centered at the origin, which means that the line is attached to the wall for —0.5 < x. A 
convergence exponent of 1.5 is observed for this formally second-order method. 


4.2.2. Spurious Dynamics of TVD Schemes for the Embid et al. Problem* 

It has long been observed that the occurrence of residual plateauing is common when TVD 
and ENO types of schemes are used to time march to the steady state. That is, the initial 
decrease in the residual levels out and never reaches the convergence tolerance. See Yee (1986, 
1989) and Yee et al. (1990) for some discussion. This has often been overcome by ad hoc 
modification of the flux limiter or similar device in problem regions. 

A recent study (Burton & Sweby 1995) investigated this phenomenon using a dynamical 
systems approach for the one-dimensional scalar test problem of Embid et al. (1984) 

u t = $(*)« = (6* — 3)u, *€(0,1), (4.9a) 

*We would like to thank Paul Burton for the computations used in this section. 
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with boundary conditions 


u(0) = 1, u(l) = -0.1. (4.9b) 

This equation with the flux function /(u) = u 2 /2 has the property that there are two entropy 
satisfying steady solutions consisting of stationary shocks jumping between the two solution 
branches 


u\x) = Zx(x — 1) + 1 (4.10a) 

u r (z) = Zx(x - 1) - 0.1. (4.10b) 

For this problem the two possible solutions consist of a single shock, either approximately at 
*i = 0.18 or = 0.82. It can then be shown (see Embid et al. 1984 for details) that the 
solution with a shock at *i is stable to perturbations while the solution with a shock at x 2 is 
unstable. 

Embid et al. solved (4.9) using three different methods — the first-order implicit upwind 
scheme of Engquist and Osher, its second-order counterpart, and the second-order explicit 
MacCormack scheme. All three schemes used time stepping as a relaxation technique for 
solving the steady-state equation. The initial conditions were taken to follow the solution 
branches (4.10) from the boundary values, with a single jump between the two branches. The 
results obtained showed that, although the implicit schemes allowed large time steps and hence 
fast convergence, if the initial jump was taken too near the unstable shock position x 2 , then for 
some ranges of Courant number. 


At 



(4.11) 


the schemes would converge to a physically unstable shock. This phenomenon was studied 
both for these three schemes and a variety of flux limited TVD schemes (Sweby 1984) in Burton 
& Sweby (1995), where not only the full problem was studied but also a reduced 2x2 system 
was investigated using a dynamical system approach. We summarize this investigation here. 


The schemes investigated were explicit and implicit versions of the Engquist-Osher and 
TVD flux limiter schemes using the minmod, van Leer, van Albada and superbee flux limiters. 
For the time discretization, forward Euler was used for the explicit implementations while 
linearized implicit Euler was used for the implicit computations. For the second-order flux 
limiter schemes the Jacobian matrix used was taken to be that of the first-order Engquist-Osher 
in order to allow easy inversion. 


The schemes are 
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(a) Explicit Scheme 

u " +1 = u" - AA -(/r +1 + /+) + Aij(«)u? 

- i A A _ [*(r+ )( A/ j+ 1 ) + - ^(r- + , )( A/,+ , )‘ ] , 

(b) Linearized Implicit Scheme 


J(“")[“? +1 - “"1 = - AA -(/ r +1 +//) + A M«K 

- iAA-[^(rt)(A/ ;+ j) + — ^( r j+, )(A/j+ J ) ] > 

where ff are the Engquist-Osher numerical fluxes 


/+ = /(max(«j, 0)), 

fj = /(min^^O)). 

The flux differences are given by 

(A/, + i)+ = -(/t +1 - f(uj+i)), (A f H ,)~ = (f; 
and the solution monitors by 


- /(“i))> 




(A/ y .)*!=“ 


,(A/ i+ t)±J 

Finally, J is the Jacobian matrix and the flux limiter <f>(r) is one of 


<^ 0 (r) = 0 — first-order Engquist-Osher (E-O) scheme, 

<f>i(r) = max(0, min(r, 1)) — the minmod limiter, 

<^ 2 (r) = max(0, min(2r ) 1), min(r, 2)) — Roe’s superbee limiter, 


<f>VL(r) 

4>VA(r) 


M- — van Leer’s limiter, 

1 + |r| 

r-fr 2 

van Albada’s limiter. 

1 + r 2 


(4.12) 


(4.13) 

(4.14s) 

(4.14b) 

(4.15) 

(4.16) 

(4.17) 

(4.18) 

(4.19) 

(4.20) 


(4.21) 
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The experiments reported in Burton & Sweby (1995) used a grid spacing of A* = 0.025 
with initial conditions consisting of a single jump between the solution branches (4.10) near 
either the stable shock (*i) or the unstable shock (* 2 )- The convergence criterion used was the 
following bound on the residual 


X>i + 1 -“il< 1 <r ,5 > (4.22) 


with an upper limit of 2000 iterations being performed. 

The results of applying these schemes to'problem (4.9) largely echoed those reported by 
Embid et al. For the explicit schemes convergence, when it occurred, was to the stable shock. 
It was found that there were regions of Courant number (c < cj) for which the schemes 
converged, regions (c< < c < Cj) for which convergence did not take place using (4.22) within 
2000 iterations, and regions (c > cj ) for which the schemes were unstable. This is summarized 
in Table 4.1. The absence of an entry in Table 4.1 corresponds to residual plateauing. 

Notice that for the superbee flux limiter there was no range of Courant numbers for which the 
scheme converged. Closer inspection reveals residual (defined as r n = |u n+1 — u n |) plateauing 
at around 10 _s . For the other schemes when Ci < c < cj the nonconvergence observed arises 
from a similar process, except that the residual does not necessarily level out completely, but 
decreases at a very gradual rate, resulting in very slow convergence. 

The implicit scheme experiments revealed that the choice of initial conditions could cause 
convergence to the unstable shock for certain ranges of Courant number. For an initial jump 
near the stable shock, the schemes (with the exception of the van Leer limiter) converge to 
the stable shock for c < 11. However, for an initial discontinuity near the unstable shock, 
convergence could sometimes be towards the unstable shock. The situation is summarized in 
Table 4.2, where again the absence of an entry corresponds to residual plateauing. 

To gain further insight into this problem. Burton & Sweby (1995) considered a reduced 
problem consisting of two free points at one of the shocks, with exact solution values being 
imposed as boundary conditions on either side. This then leads to a two dimensional dynamical 
system which, although obviously a gross simplification of the full problem, was hoped to still 
maintain some of the qualitative behavior. 

The situation is as shown in Fig. 4.6, where the free points are X and Y, the remaining 
points (Uu, Pi, Urr and U r ) being set at exact analytic values to provide boundary conditions. 
Two such values are needed on either side to provide the necessary information for the flux 
limiters. Substitution of these points into the numerical scheme then leads to a two-dimensional 
system. For example, the explicit Engquist-Osher scheme yields 
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x n+1 = X n - 


At 

40 


yn-f 1 yn 

40 


f-(Y n ) + f+(x n ) - f-(x n ) - / + (tr,) 

f-(Ur) + f + (Y n ) - f-(Y n ) - f + (X n ) 


20 


36 

20 


A tY n , 


(4.23a) 

(4.23b) 


where a step size of A* = A has been used. 

For the first-order explicit and implicit schemes some analysis on the reduced problem can be 
performed. Table 4.3 summarizes the findings. Note that for both schemes spurious fixed points 
are introduced by the simplification of the problem.' These both have X and Y of the same 
sign and would not be tolerated for the full problem. It is only the proximity of the boundary 
conditions for the reduced problem which allow them to exist as fixed points. However, the 
remaining fixed points and their stability agree well with numerical results obtained for the full 
problem. 

Analytical results could only be obtained for the first-order scheme and so numerical 
experiments were performed for the flux limiter schemes. These consisted of generating 
bifurcation diagrams for X and Y against At and the plotting of basins of attraction in the 
(X,Y) plane for fixed values of A*. The explicit schemes were shown to possess no spurious 
dynamics below their respective stability limits, apart from that introduced by the simplification 
of the problem (i.e. outside of the quadrant (X > 0, Y < 0)). As At was increased above the 
stability limit the schemes entered a period of bifurcation and chaos accompanied by a dramatic 
shrinkage in the numerical basins of attraction. 

The dynamics of the implicit schemes at the unstable shock showed the falsely stable fixed 
point becoming stable for large values of At. For all the limiters tested the stabilizing of 
the fixed point was accompanied by the introduction of additional, spurious (period 2) fixed 
points. These spurious solutions caused a reduction in size of the basin of the falsely stable 
fixed point. The fact that the more compressive limiters took longer to recover from the effects 
of the spurious fixed points seems a possible cause of the phenomenon of residual plateauing 
experienced in the full problem. Due to a space limitation, see Burton & Sweby (1995) for the 
illustrations. 

It must be realized that although the residual plateauing illustrated is around the physically 
unstable shock (to which we would usually not wish to converge), the fact that it is not a repelling 
phenomenon will in itself have repercussions on convergence to the correct, physically stable 
shock. Indeed no such nonconvergence was observed for Courant numbers greater than 1 1. We 
conclude this section by emphasizing that the reduced problem indicated a possible cause for 
residual plateauing, although it should be realized that the dynamics of the full problem does 
not coincide precisely with that of the reduced problem. 
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4.2.3. The Dynamics of Grid Adaptation 

Consider a model convection-diffusion equation of the form 

ti t + /(«)• = «*•» 


(4.24) 


with the linear case, /(«) = u and the nonlinear case, /(«) = f u 2 (the Burgers equation). The 
boundary conditions for the linear case are u(0,f) = 0 and u(l,t) = 1 and for the nonlinear 
case, 1 and u(l, t) = -1. For small e, these boundary conditions result in steady-state 

solutions of a boundary layer at * = 1 and a viscous shock at * = 1/2, respectively. In both 
cases the steepness of the feature is governed by the diffusion coefficient parameter «. Besides 
its steepness feature, one of the main reasons for considering the linear convection-diffusion 
equation is to show that grid adaptation alone and/or nonlinear schemes such as TVD schemes 
can introduce unwanted dynamics to the overall solution procedure. The authors realize that 
model (4.24) is not the best model to illustrate the dynamics of the studied schemes since the 
model is not stable under perturbations. However, it serves to show what type of spurious 
numerics would occur under such an environment. 

One common criterion used for grid adaptation is the equidistribution of a positive definite 
weight function to(z,t), often taken to be some monitor of the numerical solution u(z,t) of 
the underlying PDE. A grid zo < zj (<)<••• < zj_i(i) < xj, where zo and zj are fixed, 
equidistributes tt>(z, i) (at time t) if 

J w(x,t)dx = J w(x,t)dx = j j w(x,t)dx , (4-25) 

for j = 1, . . . , J. A one-parameter family of weight functions 


w(x,t) = -y/l - a + cru*(z,t), a £ [0,1], 


(4.26) 


can be chosen where a — 1/2 corresponds to equidistribution in the arc-length and a = 0 
yields a uniform grid. Approximating w(x, t) to be constant in each interval (z,-_i, xj) yields 


(4.27) 


(-y/l - a + at {xj - Xj- 1) = (-y / 1 - a + au «) J+ i/j ( Xj+1 ~ 

Given a numerical solution of the PDE we can approximate the derivatives by 

«.| j-i,2 * (4.28) 

Xj- Xj- X 

Equation (4.27) is nonlinear in {z ; } if we use (4.28). However, (4.27) is linear in {z,} if {xj} 
in (4.28) uses the existing grid. In this case we can solve the tridiagonal system (4.27) for a new 
set of {xj} to obtain an updated grid. 
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Given a set of initial data and an initial grid, the procedure is to numerically solve the PDE 
and (4.27) in a time-lagged manner. We use nodal placement and the l 2 norm of the solution 
to illustrate our results. We use the previous time step value for »j in (4.28) to achieve a linear 
tridiagonal system for the updated grid in (4.27). Our preliminary study shows that the solution 
procedure of Ren and Russel (1992), and Budd et al. (1995a) in solving the coupled PDE and 
(4.27) are less stable than the present linearized form. See also Neil (1994) for a similar study 
and conclusion. The regridding strategy adopted here was to regrid after every time step of the 
PDE method, either interpolating updated solution values from the old grid or performing no 
adjustment at all due to grid movement. This latter approach in effect presents the PDE method 
with new initial data to the problem at each step. 

The dynamics of the above one-parameter family of mesh equidistribution schemes coupled 
with different spatial and time discretization were studied numerically in Sweby & Yee (1994) 
and Yee & Sweby (1995a,b) using the above numerical procedure. The spatial discretizations 
included three-point central, second-order upwind and second-order TVD schemes. The time 
discretizations included explicit Euler, second- and fourth-order Runge-Kutta and the linearized 
implicit Euler methods. In a parallel study, Budd et al. (1995b) made use of the AUTO computer 
bifurcation package (Doedel 1986) to obtain bifurcation diagrams for similar grid adaptation 
methods for the steady part of the above PDEs with a different form than (4.27). However, 
the dependence on known solutions of the discretized PDEs and grid equations as starting 
values limits its usage. In Sweby & Yee (1994) we utilized the power of the highly parallel 
Connection Machine CM-5 to undertake a purely numerical investigation into the dynamics of 
the time-marching adaptive procedure. 

We divided a chosen parameter space (e.g., e) into 512 equal intervals, with all other 
parameters (a, At, initial data) fixed. For each chosen parameter value, we iterated the 
discretized PDE and the grid function, in general, 4000 steps (8,000 steps for explicit methods) 
to allow the solution to settle to an asymptotic state. Then we performed a series of time 
step/regridding stages, during which we investigated the dynamics by producing an overlaid 
plot of the l 2 norms of the numerical solution and the grid distribution at each step. This 
resulted in a bifurcation type diagram or the grid displacement diagram as a function of the 
physical or discretized parameters. We also performed numerical studies by only preiterating 
the discretized PDEs to the steady state for a fixed grid before solving both the discretized PDE 
and the grid adaptation function. We found in most cases this solution process is less stable and 
more likely to get trapped in a spurious mode than in the fully coupled process. 

For this study, we took into consideration the grid density, an even and odd number of nodes, 
and whether or not there is interpolation after each regridding. The grid density studies consist 
of 4, 5, 6, 9, 10, 19, 20, 49 and 50 grid points. There is no apparent sign of even or odd grid 
dependence. The resolution and stability of TVD schemes are also grid independent. However, 
the central difference scheme experienced instability more often for coarser grids, and the 
second-order upwind is slightly more stable, with better resolution than the central scheme. As 
expected, the stable time step required for the explicit methods was orders of magnitude lower 
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than that for the implicit method. For the TVD schemes, comparison of the dynamical behavior 
of the five limiters of (3.50) of Yee (1989) was performed. Four of the limiters are the same as 
(4.18)-(4.21). Due to the simplicity of the PDEs, their dynamical behavior is similar, although 
there were slight differences in the stability and resolution. Due to space limitations, we 
summarize the results without presenting the actual computational figures. Interested readers 
should refer to our original papers for details. 

Summary : We consider separately the cases “with” and “without” interpolations. The term 
“scheme” from here on means the overall adaptive scheme procedure. 

(a) No Interpolation: For the linear problem, the behavior of the adaptive TVD schemes is 
similar to that of the classical shock-capturing methods. As opposed to the uniform grid case, 
the adaptive TVD schemes without interpolations behave rather poorly in terms of stability and 
allowable e values. See Figs. 1 1 and 12 of Yee & Sweby (1995b). The solutions refuse to settle 
down for larger A t and/or smaller e. For the nonlinear problem, the behavior of the adaptive 
TVD schemes is similar to that of the uniform grid case. The range of allowable e and At in 
terms of stability and convergence rate and settling of the grid distribution are far better than in 
the linear problem. It appears that for problems with shocks, adaptive TVD schemes prefer no 
interpolation after each regridding. See Figs. 13 and 14 of Yee & Sweby (1995b). 

(b) With Interpolations: For the linear problem, as expected, both adaptive TVD schemes 
and adaptive classical schemes behave in a similar manner in terms of stability and convergence. 
Adaptive TVD schemes are less stable and have a smaller allowable range of e than the uniform 
grid case. Overall, adaptive TVD schemes behave far better than their counterparts without 
interpolation for the linear problem as can be seen in Figs. 1 1 and 12 of Yee & Sweby (1995b). 
For the nonlinear problem, the adaptive TVD schemes with interpolations behave like the 
classical shock-capturing method. They experience nonconvergence of the solution, and the 
grid distribution cannot settle down for a certain range of At. This can be seen in Figs. 13 and 
Hof Yee & Sweby (1995b). 

It is surprising to see the opposite behavior of the adaptive implicit TVD schemes for the two 
model PDEs with and without interpolation combinations, especially when the same physical 
parameters, discretized parameters and initial data were used. 


4.3. Mismatch in Implicit Schemes for Time-Marching Approaches 

When implicit methods are used to time-march to the steady states, it is a common practice 
in CFD to use a linearized and/or a simplified implicit operator (or mismatched implicit/explicit 
operators) to reduce operations count. The simplified implicit operators usually are first-order 
and the explicit operator retains higher order spatial accuracy. The simplified implicit operator 
might not be consistent with the original implicit scheme. It might also be nonconservative even 
though the original and/or the explicit operators are conservative. One popular formulation 
with these mismatched implicit/explicit operators is the “delta formulation” of Beam & 
Warming (1978) in conjunction with implicit LMMs time discretizations. The logic is that 
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if the solution converges, the explicit operator dictates the final accuracy of the steady-state 
numerical solution. As discussed in Section 2.4, these mismatched implicit/explicit operators 
might induce unwanted spurious dynamics into and/or reduce the convergence rate of the 
overall solution procedure. To illustrate just this point, we summarize some old work of Mulder 
& van Leer (1984) and the first author’s experiences (Yee 1986, 1989, 1990) in selecting the 
more desirable mismatch operators. These works use the delta formulation with a variant of 
the first-order upwind spatial discretizations for the implicit operator. The time discretization is 
the implicit Euler with the noniterative linearized form as discussed in Section 3.5.2. After the 
linearization in time and the drop in spatial discretization to first order, there are many ways to 
approximately evaluate the Jacobian matrix associated with the linearization for high-resolution 
and higher-order upwind shock-capturing schemes. Due to a space limitation, the readers are 
referred to the original papers for details or Yee (1989) for a summary. 

Mulder and van Leer studied two first-order upwind spatial discretizations (van Leer’s 
differentiable and Roe’s nondifferential forms) for the implicit operator and a first-order or 
second-order upwind spatial discretization for the explicit operator. They concluded that the 
differentiable first-order upwind implicit operator gives quadratic convergence for flow of an 
isothermal gas along an almost circular path through the stellar gravitational field of a rotating 
two-armed spiral galaxy. The governing equations are a system of hyperbolic conservation 
laws. With Roe’s nondifferentiable split flux-differences the iterations may get trapped in a 
limit-cycle. For a comparison of the two implicit operators, they used the same grid, physical 
parameters, explicit operator, initial data and numerical boundary treatment. 

Yee (1986) constructed conservative and nonconservative linearized forms of implicit TVD 
schemes. A rather detailed study on the convergence properties of these forms was performed. 
Both the linearized nonconservative implicit (LNI) and the linearized conservative implicit 
(LCI) forms were of first or second-order accurate in space. Numerical experiments in Yee et 
al. (1985), Yee (1986), Yee & Harten (1987) and Yee (1989) showed that both the second-order 
LNI and LCI forms are very unstable even if a very small time step is used. The first-order 
LNI and LCI perform far better. However, the first-order LCI form is more stable and has 
a better convergence rate than the LNI counterpart. The first-order LNI form has a higher 
chance of converging to a nonphysical solution. Also the residual sometimes stagnates or gets 
trapped in a limit cycle. These conclusions are based on comparing the two forms for a variety 
of one-dimensional and two-dimensional practical problems containing complex shock waves. 
For both LCI and LNI forms, the more compressive limiters, such as the superbee, are very 
unstable. The residual stagnates even with very small time steps. See Yee et al. (1990) for 
their performance for hypersonic computations. Grid refinement in this case does not improve 
the situation. Again, in order to isolate the cause of the convergence problem, the same grid, 
physical parameters, explicit operator, initial data and numerical boundary condition treatment 
were used. In passing, both the first-order LCI and LNI are heavily used in the CFD community 
with TVD schemes other than the Harten & Yee type (Yee 1989). This includes but is not 
limited to the various UNO, ENO and high-order upwind schemes for the explicit operator. 

The above convergence phenomena using pseudo time marching approaches in conjunction 
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with high-resolution implicit TVD schemes contributed to yet another kind of spurious numerics, 
which are very different from that in Sections 4. 1 and 4.2. From a dynamical systems standpoint, 
they illustrate the fact that before a steady state is reached, the nonlinear difference equation 
representing each of these simplified implicit operators are distinct discrete dynamical systems. 
Before a steady state has been reached, during transient states, the solution procedures take 
different paths to get to the steady state, depending on the implicit operator, the time steps, 
initial conditions, and grid spacings (even with the same explicit operator). Some combination 
of these choices can get trapped in a stable spurious numerical solution. Other combinations 
of these choices can by pass these traps. That is why there are so many forms of the implicit 
operator being used. Some implicit operator perform better than others, depending on the 
physical problem. Even after a steady state has been reached and the residual error of the 
explicit operator is zero, the solution can still be spurious, since a stable spurious steady state 
would produce a machine zero residual error if the spurious behavior is due to the spatial 
discretization (see Section II under the sub-heading “Reliability of Residual Tests”. 


V. Spurious Dynamics in Unsteady Computations 

In Section 3.8, we cited that numerics can not only introduce and suppress chaos but 
also introduce chaotic transients. Here, some examples from CFD computations that exhibit 
analogous spurious behavior are illustrated. Sections 5.1, 5.2, 5.5 and 5.6 are written by the 
authors of the original work. Shi Jin of the Georgia Institute of Technology summarizes his 
recent work on oscillations induced by numerical viscosities in Section 5.1. Shi’s work has 
an important implication in spurious dynamics using time-marching approaches as well as 
slowly moving shock waves in transient computations. Section 5.2. summarizes the work 
of Brown & Minion (1995). It is concerned with spurious vortices in two-dimensional thin 
shear layer incompressible flow simulations using high-resolution shock-capturing methods. 
Laurence Keefe of Nielsen Engineering summarizes his unpublished work on chaotic transient 
computation that he performed in the late 1980’s in Section 5.3. Sections 5.4 - 5.6 are the joint 
work of the first author with John Torczynski of Sandia National Laboratories, Scott Morton and 
Miguel Visbal on spurious behavior of underresolved grid and semi-implicit time-discretizations 
(Yee et al. 1997). For additional results concerning other issues, see Moore et al. (1990), 
Corless (1992, 1994), Poliashenko & Aidun (1995) and Read & Thomas (1995). 

5.1. Oscillations Induced by Numerical Viscosities in 1-D Euler Computations 


Earlier work has reported the difficulty of computing slowly moving shocks (Robert 1990, 
Woodward & Colella 1984), where first-order Godunov (Gounov 1959) or Roe-type methods 
produce spurious long wave oscillations behind the shock and eventually contaminate the 
downstream pattern. Here slowly moving means that the ratio of the shock speed to the 
maximum wave speed in the domain is much less than one. Several heuristic arguments, or 
improvements on the Riemann solver have been made in Arora & Roe (1996), Jin & Liu 
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(1996), Liu (1995) and Woodward & Colella (1984). To investigate the dynamical behavior of 
shock-capturing methods for slowly moving shocks, we review the work of Jin and Liu (1996) 
using the traveling wave analysis and stability theory of discrete shocks. The goal is to carefully 
study this peculiar numerical phenomenon and to understand its formation and propagation, 
instead of solving this problem. 

Recall that the definition of a discrete traveling wave solution an approximation of 
U(xj,t n ),t n = nAt, requires 

, (5.1) 

3 J-np’ v ’ 

where 

iAt/Ax = p/q, (5-2) 


with * the wave speed for some relative prime integers p and q. During a numerical calculation 
condition (5.2) may not hold. In other words, at different times the numerical viscous profiles 
correspond to different families of the traveling waves. Downstream oscillations are generated 
by such perturbations of the discrete traveling wave profile. The oscillations propagate 
along characteristics and behave diffusively (decay in Li and Loo). Th e perturbing nature 
of the viscous shock (traveling wave) profile is the constant source for the generation of the 
downstream oscillations for all time. In their numerical experiments, Jin & Liu also observed 
the periodic structure of the perturbing viscous shock profile. The period is essentially the time 
for the shock to propagate one spatial grid cell. 

In Jin & Liu (1996) a numerical example, using a Roe-type upwind scheme on a Riemann 
problem of the compressible Euler equations that admits slow shocks was given. Among the 
numerical artifacts observed in that example are the momentum spikes at the shock, and the 
downstream oscillations. They are indeed numerical artifacts. The momentum spike is generated 
by the artificial numerical viscosity introduced in the continuity equation, which does not exist 
in the physical Navier-Stokes equations. The downstream oscillations are introduced by the 
dynamical behavior of the numerical viscous traveling wave profile. To relate this phenomena 
with dynamical systems, a traveling wave analysis on a “viscous isentropic Euler equations” 
formulation (Euler equations with a special linear viscosity term in both the continuity and 
momentum equations) is presented in Section 5.1.1 to show the existence of the momentum 
spike. This is compared with the momentum profile of the Navier-Stokes equations, which does 
not have the spike. For the dynamics of the downstream oscillations and its relation with the 
stability and perturbation of the discrete shocks see Jin & Liu (1996) for details. 

5.1.1. The Momentum Spikes 

This section present a traveling wave analysis of the “viscous Euler equations” in the 
formation of the momentum spike. Consider the following “special viscous isentropic Euler 
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equations” for density p and momentum m: 

d t p + d,m = td ma p, (5.3a) 

d t m + d m h p(p) = ed Ba m. (5.3b) 

. P 

Here the pressure p(p) = kp y for some constants k and 7 . This hyper bolic sys tem has two 
distinct eigenvalues u ± c, where u = m/p is the velocity and c = Vr is the sound 
speed. Although the true numerical viscosity is far more complicated than those appearing 
on the right-hand-side of (5.3), a study of (5.3) is sufficient for a full understanding of the 
numerical momentum spike. 

The traveling wave solution to (5.3) is examined. Let £ = * e tt , where s is the shock speed. 
Then the traveling wave solution takes the form 

p(z,t) = <!>(£ ), m(aj,t) = ij>(£) (5.4) 



(5.9b) 
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This system has two fixed points: V+ = (<p+,tp~) on the right and V- = (</>-,*/>-) on the 
left in the phase plane of (<p,ip). It has two distinct eigenvalues, Ai = u - c and A 2 = u + c, 
with corresponding eigenvectors Ri = (l,t* — c) T and Rj = (l,t* + c) T . By the entropy 
condition (5.8), V+ is a saddle point with a stable manifold on Ri, and V- is a source. Thus a 
heteroclinic orbit O will connect V- and V+ in the direction of Ri (Smoller 1983), as shown in 
Fig. 5.1. (In Fig. 5.1 R* and Rf are the two eigenvectors at V± respectively). The orbit O is 
smooth, and ^ is not identically zero if <p~ ^ . Thus (5.9a) implies that tp is not a constant. 

Moreover, whenever <p(() connects <p~ and <p+ with a monotone profile, <p becomes a spike. 
Thus ip = j> + ip- is a spike. 

For a nonstationary shock, the traveling wave solution (5.4) applied to (5.3a) gives 

4> = —s(<p — <p.~) - . (5.10) 

or 

ip = s<p + ^ + (t/>_ — (5-11) 

Hence tj> is a superposition of a monotone profile s<f> with a spike corresponding to 4>. When s 
is small (for a stationary or a slowly moving shock), the monotone profile s<f> becomes small 
and the spike term <f> dominates. Thus the shock profile of i/> is a non-monotone spike. Therefore 
the spike is usually generated in a stationary or slowly moving shock, as shown in the earlier 
examples. For a strong shock the monotone profile s<f> dominates so the shock profile of the 
momentum is monotone. 

Since the more physical viscous shock profile is determined by that of the Navier-Stokes 
equations, the viscous profile of the isentropic Navier-Stokes equations is now studied and 
compared it with that of the viscous Euler equations (5.3). The isentropic Navier-Stokes 
equations are 

dtp + d u m = 0, (5.12a) 

d t m + d 9 h p(p) = cd„ m • (5.12b) 

Applying the traveling wave solution (5.4) in (5.12) and again assuming the shock speed s 0, 
one obtains the following ODEs: 

ip = t/>_ , (5.13a) 

(^) = ^ (5.13b) 

Equation (5.13a) shows that tp is a constant and thus contains no spike. Let f(<p) = -j- 
then (5.13b) becomes 

(|) =m -w-)- 


(5.14) 
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Since f(<f>) is a strict convex function and = /(^+), f{<f>) ~ is always negative 

between <p- and <p+. Therefore, ^ does not change sign by (5.14). This implies the 

monotonicity of l/<f> or <p. 

When s ^ 0, applying the traveling wave solution (5.4) to (5.12a) gives 

ip = s<p + {ip- — *<!>-)’ (5.15) 

Thus whenever <p is monotone so is ip. This excludes the possibility of a momentum spike for a 
moving viscous shock in the Navier- Stokes equations. 

In conclusion, even if the viscous profile of <p of the Navier-Stokes equations (5.12) could 
be similar to that of the viscous Euler equations (5.3), the profiles of ip may be significantly 
different. Since the physically relevant solution of the Euler equations is considered to be 
the zero viscosity limit of the Navier-Stokes equations, the momentum spike appearing in the 
viscous Euler equations is totally nonphysical. 

By examining the interrelation between the viscous Euler and the Navier-Stokes equations 
one can come up with a change of variables that recovers the Navier-Stokes equations in 
an asymptotic sense from the Euler equations. This also motivates a numerical change of 
variable from the cell-center or cell average momentum to the mass flux, which eliminates the 
momentum spike exactly. For details see (Jin & Liu 1996). However, what is more catastrophic 
is the downstream oscillations, which cannot be easily eliminated. See Jin & Liu (1996) for 
details. 


5.1.2. Discussions 

Behaviors similar to that studied in Jin & Liu (1996) occurs in schemes that are of monotone, 
TVD or ENO type. Note that all these monotonicity theories are established only for scalar 
equations, or linear systems via the characteristic variables. For nonlinear systems there are no 
global characteristic variables. Thus these methods are usually extended to nonlinear systems 
using the idea for linear systems; i.e., via the so-called local characteristic decomposition (using 
the Roe matrix (Roe 1981) example). Since there is no theory for the monotonicity of these 
methods for nonlinear systems, it is not surprising to see the non-monotone behavior represented 
by the spike and the downstream oscillations reported here. Jin & Liu (1996) pointed out that 
to fully solve this problem, one may need a method that is systematically “monotone, TVD or 
ENO” instead of applying a scalar monotone, TVD or ENO scheme to nonlinear systems. One 
may also need to choose numerical viscosity properly so it mimics the physical viscosity of the 
Navier-Stokes equations. The ultimate goal is to have a scheme that not only provides a high 
resolution but, more importantly, has a more stable viscosity profile. 

The novelty of the work of Jin and Liu is that they are the first to use the traveling wave 
analysis to prove the non-monotonicity of the solution for nonlinear systems, and to link the 
downstream oscillations to the stability of discrete traveling wave profiles. To really understand 
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the behavior of shock-capturing methods for nonlinear systems, and to ultimately design 
nonoscillatory schemes for nonlinear systems, good theories for both inviscid and viscous 
nonlinear systems need to be developed. This remains an open and challenging research subject 
for the future. 

5.2. Spurious Vortices in Under-Resolved Incompressible Thin Shear Layer Flow Simulations 

Brown & Minion (1995) performed a thorough study of a second-order Godunov-projection 
method and a fourth-order central difference method for the two-dimensional incompressible 
Navier-Stokes equations as a function of the resolution of the computational mesh, with the rest 
of the physical and discretized parameters fixed. In the authors’ opinion, this is a good example 
of isolating the cause of the spurious behavior. The physical problem is a doubly periodic 
double shear layer. The shear layers are perturbed slightly at the initial time, which causes the 
shear layer to roll up in time into large vortical structures. For a chosen shear layer width that is 
considered to be thin and a fixed perturbation size, they compared the solution for four different 
grid sizes (64 x 64, 128 x 128, 256 x 256, 512 x 512) with a reference solution using a grid 
size of 1024 x 1024. For the 256 x 256 grid, a spurious vortex was formed midway between 
the periodically repeating main vortex on each shear layer. The 128 x 128 solution showed 
three spurious vortices along the shear layer. The spurious vortex disappeared with a 512 x 512 
mesh. They also disabled the flux limiters (a strictly upwind Fromm’s method), and found 
the behavior to be similar. A subsequent study (Minion & Brown 1997) using five different 
formulations and six different commonly used schemes in CFD found similar behavior. They 
concluded that the spurious vortex is the artifact of underresolution of the grid and the behavior 
is caused by a nonlinear effect.. Linking this behavior with a re-interpretation of their conclusion 
using nonlinear dynamics, we interpret their observation as follows. For the particular grid size 
and time step combination, stable spurious equilibrium points were introduced by the numerics 
into a portion of the flow field while the major portion of the flow field was predicted correctly. 

In other words, the spurious vortices are the solution of the discretized counterpart for that 
particular range of grid size and time step. The number of stable spurious vortices is a function 
of the grid size. As the grid spacing decreases, the spurious equilibriums gradually become 
unstable and the numerical solution mimics the true solution. 

5.3. Chaotic Transients Near the Onset of Turbulence 
in Direct Numerical Simulations of Channel Flow 

5.3.1. Motivation 

In addition to the inherent chaotic and chaotic transient behavior in some physical systems, 
numerics can independently introduce and suppress chaos as well as chaotic transients. A brief 
background on chaotic transients and the significance and implications of Keefe’s work on 
numerical uncertainties is needed before presenting his results. Loosely speaking, a chaotic 
transient behaves like a chaotic solution (Grebogi et al. 1983). A chaotic transient can occur in 
a continuum or a discrete dynamical system. One of the major characteristics of a numerically 
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induced chaotic transient is that if one does not integrate the discretized equations long enough, 
the numerical solution has all the characteristics of a chaotic solution. The required number 
of integration steps might be far beyond those found in the standard CFD simulation practices 
before the numerical solution can get out of the chaotic transient mode. Furthermore, standard 
numerical methods, depending on the initial data, usually experience drastic reductions in step 
size and convergence rate near a bifurcation point of the continuum in addition to the bifurcation 
points due solely to the discretized parameters. See Yee & Sweby (1996a, 1994, 1995a) for 
a discussion. Consequently, the possible numerically induced chaotic transient is especially 
worrisome in direct numerical simulations of transition from laminar to turbulent flows. Except 
for special situations, it is impossible to compute the exact transition point by mere DNS of the 
Navier-Stokes equations. Even away from the transition point, this type of numerical simulation 
is already very CPU intensive and the convergence rate is usually rather slow. Due to limited 
computer resources, the numerical simulation can result in chaotic transients indistinguishable 
from sustained turbulence, yielding a spurious picture of the flow for a given Reynolds number. 
Consequently, it casts some doubt on the reliability of numerically predicted transition points 
and chaotic flows. It also influences the true connection between chaos and turbulence. The 
present section makes use of this knowledge from continuum and discrete dynamical systems 
theory to identify some of the aforementioned numerical uncertainties. 

5.3.2. Numerical Simulation 

Numerical simulations of wall-bounded turbulent shear flows have proven to be a powerful 
tool for investigating both the physics and mathematics underpinning these practically important 
flows. Boundary layers and turbulence are inescapable in aeronautical applications, and simple 
flows over a flat plate or within a channel provide an accessible arena in which to understand 
turbulent dynamics and to test ideas for its modification and control that can benefit the 
performance of real aircraft. 

Direct numerical simulations (DNS) of turbulent channel flow are the best developed of 
these techniques, with a better than 20 year history (Deardorff 1970, Schumann 1973, Moin & 
Kim 1982), and a relatively quick maturity (Kim et al. 1987) due to the favorable mapping of 
high accuracy spectral methods onto the geometry and known phenomenology of this flow. The 
physical situation is depicted in Fig. 5.2, where a flow is confined between planes at y = ±1 
and is driven in the * -direction by a mean pressure gradient dp/dx. The flow is characterized 
by a Reynolds number Re = U^L/v, where Uoo is the mean centerline velocity, L is the 
channel half-height, and v is the kinematic viscosity. Within the channel the flow satisfies 
the incompressible Navier-Stokes equations and the no-slip boundary conditions are applied 
at the walls. In the particular calculations shown here these equations have been manipulated 
into velocity- vorticity form, where one integrates equations for the wall-normal velocity v, and 
normal vorticity rj and recovers the other two velocity components from the incompressibility 
condition and definition of rj. 


£a‘v = fc. + 


(5.16a) 



71 


where 


0 

at 


dx 


' = *• + jb A2 ’ 

(5.16b) 

O 

II 

Cfc |C& 

+ 

(5.16c) 

dw du dw 

dz y ^ dz dx 

(5.16d) 

+ _ 8r) + («? + a^) ffj 

(5.16e) 

- dHi 9Hi 
8z dx' 

(5.16f) 


Here the Hi contain the nonlinear terms in the primitive form of the Navier-Stokes equations 
and the mean pressure gradient. 


Two experimentally observed facts figure strongly in the choice of numerical method used 
to integrate the Navier-Stokes equations in this geometry: the velocity increases extremely 
rapidly normal to the wall, and turbulent channel flows are essentially homogeneous in planes 
parallel to the wall. The first requires a concentration of grid points near the wall, and the 
second suggests use of a doubly periodic domain in planes parallel to the wall. Happily, a high 
accuracy spectral representation of the velocity field (u, v,to) meets these needs: 


s = E £ £ (5-17) 

2 m n 

where the Ti(y) are Chebyshev polynomials. The numerical problem then becomes dependent 
on a and /3 in addition to Re. For the time discretization, mixed explicit-implicit methods are 
used. The nonlinear terms in the equations are advanced using second-order Adams-Bashforth 
or a low storage, third-order Runge-Kutta scheme (Spalart et al. 1991), while the viscous terms 
are advanced by Crank-Nicholson. The algorithm based on this spectral spatial representation 
and mixed time advance has been tested and used extensively, demonstrating an ability not 
only to reproduce experimental results, but to go beyond them in elucidating flow features not 
easily investigated in experiments. This code, and ones similar, (Handler et al. 1989, Jung et al. 
1992) are currently being used to investigate a variety of turbulence control ideas suitable for 
turbulent drag reduction and (conceivably) separation control. 

One of the central problems in studies of wall bounded shear flows is the determination 
of when a steady laminar flow becomes unstable and transitions to turbulence. In dynamical 
systems’ terms, the Navier-Stokes equations always have a fixed point solution for low enough 
Reynolds numbers, but for each flow geometry the Reynolds number at which this fixed point 
bifurcates needs to be determined. In channel flow the fixed point solution (a parabolic velocity 
profile across the channel, u(y) = (1 - y 2 )) becomes linearly unstable at Re = 5, 772 (Orszag 
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1971). However, since turbulence appears in experiments at much lower Reynolds numbers, 
it was conjectured that this bifurcation must be subcritical. Subsequent numerical solution 
of the nonlinear stability equations (Herbert 1976, Ehrenstein & Koch 1991) demonstrated 
this to be true, showing that limit cycle solutions with amplitude e branch back to lower 
Reynolds numbers before subsequently passing through a turning point and curving back 
toward higher Reynolds numbers. Thus for Reynolds numbers just above the turning point the 
flow equations have at least four solutions: the fixed point; two unstable limit cycles; and a 
chaotic solution (experimentally observed turbulence). Determining the location of the turning 
point in (a,/3,e, Re) space is known as the minimum-critical-Reynolds-number problem, and 
its solution is by no means complete. 

One way to investigate the turning point problem is to perform DNS of channel flow for 
conditions believed to be near this critical condition. Beginning with a known turbulent initial 
condition from higher Reynolds number, one integrates the flow solver in time at the target 
Reynolds number to determine whether the flow decays back to the fixed point or sustains itself 
as turbulence. Although this may not be the most efficient way to bracket the turning point, 
it has the advantage that the peculiar dynamics of the flow near the turning point, whether in 
decay or sustained turbulence, are observable. This yields information about the path along 
which flows become turbulent at these low Reynolds numbers. 

Unfortunately the flow dynamics are very peculiar near the turning point, and extremely 
long chaotic transients are observed in the computations that make a fine determination of that 
point all but impossible by this method. This can be seen in Fig. 5.3, where a time history of 
the turbulent energy in a channel flow (energy above that in the laminar flow) is plotted for a 
Reynolds number of 2, 1 9 1 . To understand the time scale of the phenomenon some experimental 
facts need to be recalled. In typical experimental investigations of channel flow, the infinite 
transverse and streamwise extent of the ideal flow are approximated by studying flow in high 
aspect ratio (10-40) rectangular ducts that typically are 50-100 duct heights long. If times 
are non-dimensionalized by the centerline mean velocity Uoo and the duct half height L, then 
statistics on turbulence are gathered by averaging hot-wire data over intervals AtUoo/L ~ 200. 
In the simulations and figure the time scale is based on the friction velocity tt T and L, where 
typically 15-20 u T ~ Uoo- Thus averaging over intervals A tu r /L ~ 10 should and does yield 
stable flow statistics that compare well with experiments. This can be seen in Figs. 5.4-5.6, 
where the near-wall velocity profile, cross- channel turbulence intensities, and Reynolds and 
shear stress distribution for the A tu r /L ~ 10 interval near the end of the transient, delineated 
by the arrows in Fig. 5.3, are shown. In each case they correspond well to available experimental 
data. Yet look at the time scale of the transient; it spans A tu r /L ~ 300 , thirty times longer than 
the time needed to obtain stable statistics that would convince most experimentalists that they 
are viewing a fully developed turbulent channel flow. This is further complicated by the wide 
variation of the transient length, dependent upon both the grid resolution (number of modes in 
the spectral representation) and the linearly stable time step of the integration. In fact, for fixed 
(a,/?, Re) it is possible to obtain sustained turbulence for one time step, but see it rapidly decay 
to the laminar flow for another, lower value of the step. 



5.3.3. Discussions 


Extended chaotic transients near bifurcation points are not an unknown phenomenon; the 
"meta-chaos" of the Lorenz system is but one of many known examples. However, the 
practicalities of numerical computation in fluid dynamics usually interfere with ones ability 
to discern whether a transient, or sustained turbulence, is being calculated. The computations 
required to obtain the transient plot in Fig. 5.3 needed 40 hours of single processor time on 
a Cray XMP, some ten years ago. Such a small amount of expended time was only possible 
because the spatial resolution of the calculation was relatively coarse (32 x 33 x 32), in keeping 
with the large scales of the phenomena expected at these flow conditions. Higher resolution 
calculations (192 x 129 x 160) (Kim et al. 1987) at greater Reynolds numbers typically have 
taken hundreds of hours (~ 250) to barely obtain the A tu r /L = 10 averaging interval that is 
so inadequate for detecting transients. Because such calculations are so time consuming, one 
typically chooses an integration time step that is a substantial fraction of the linear stability 
limit of the algorithm, so as to maximize the calculated "flow time" for expended CPU time. 
However, it is clear from these transient results that this practice has some dangers when close 
to critical points of the underlying continuous dynamical system. Thus it appears that just as 
pseudo-time integration to obtain steady solutions can result in spurious results, genuine time 
integration can result in chaotic transients indistinguishable from sustained turbulence, also 
yielding a spurious picture of the flow for a given Reynolds number. 

To conclude this section, we give a feel for the number of integration steps required for the 
above numerical solution to get out of the chaotic transient mode. This transient calculation was 
performed using a time step of .0025. In these same units (time scaled by wall friction velocity 
u T and channel half height L ) the transient calculation length was 409.6. Thus this calculation 
extended over 163,840 time steps. 

In Keefe et al. (1992) the dynamics of the computation in terms of “dimension” and 
“Lyapunov exponent” calculations were performed at a higher Reynolds number Re = 3200. 
For this higher Reynolds number, a transient occurred using a smaller time step of 0.0015 and 
a transient calculation length of 644.52 or 429,680 time steps. To examine if chaos occurred, 
Keefe et al. then took a central portion of that calculation at around 45,600 time steps and 
computed the Lyapunov exponent hierarchy from it. See Keefe et al. for additional details. 


5.4. Temporal & Spatial Refinement Studies of 2-D Incompressible Flow Over a Backward 
Facing Step 


5.4.1. Background and Objective 

The 2-D incompressible flow over a backward facing step has been addressed by many 
authors using a wide variety of numerical methods. Figure 5.7 shows the flow geometry. 
Fluid with constant density p and viscosity p enters the upstream channel of height h with a 
prescribed velocity profile (usually parabolic). After traveling a distance 1, the fluid passes over 
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a backward-facing step of height a and enters the downstream channel of height H = h + s. 
After traveling a distance L downstream of the step, the fluid exits the region of interest. For 
Reynolds numbers considered here, the flow separates at the comer and forms a recirculating 
region of length Xr behind the step. Additional recirculating regions form on the upper and 
subsequently the lower walls of the downstream channel as the Reynolds number is increased. 

Results of sustained unsteady flow from various numerical simulations have been reported 
for Reynolds numbers (Re) ranging from 250 up to 2500. The formulations included the 
vortex method, unsteady equations in stream function form, steady equations and the associated 
linear-stability problem, and the unsteady equations in primitive variable form. The numerical 
methods used cover almost all of the existing schemes in the literature. The majority of the 
numerical results are summarized in Gresho et al. (1993). The work of Gresho et al. was an 
answer to a controversy concerning the stability of the stationary solution at Re = 800. It was 
concluded by Kaiktsis et al. (1991) that transition to turbulent flow has occurred at Re = 800. 
Kaiktsis et al. examined the long-time temporal behavior of the flow and found that the flow is 
steady at Re = 500, time-periodic at Re = 700, and chaotic at Re = 800. Gresho et al. did a 
detailed grid refinement study using four different numerical methods and concluded that the 
backward facing step at Re = 800 is a stable steady flow. 

In addition to the study of Gresho et al., an extensive grid refinement study of this flow using 
a spectral element method was conducted in Torczynski (1993). The simulated geometry and the 
numerical method corresponds to that of Kaiktsis et al. (1991). Flow was examined at Reynolds 
numbers of 500 and 800. His systematic grid refinement study was performed by varying both 
the element size and the order of the polynomial representation within the elements. For both 
Reynolds number values with the transient computations stopped at t = 800, it was observed 
that low-resolution grid cases exhibit chaotic-like temporal behavior whereas high-resolution 
grid cases evolve toward asymptotically steady flow by a monotonic decay of the transient. The 
resolution required to obtain asymptotically steady behavior is seen to increase with Reynolds 
number. These results suggest that the reported transition to sustained chaotic flow (Kaiktsis 
et al., 1991) at Reynolds numbers around 700 is an artifact of inadequate spatial resolution. 
Torczynski’ s conclusion was further confirmed by a subsequent study of Kaiktsis et al. (1996) 
and Fortin et al. (1996). Fortin et al. employed tools from dynamical systems theory to search 
for the Hopf bifurcation point (transition point). They showed that the backward-facing step 
remains a steady flow at least up to Re = 1600. The purpose of the present study is to refine 
the study of Torczynski (1993) using dynamical systems theory to interpret the results. The 
next two sections give details of Torczynski’ s 1993 analysis and the present study. 

5.4.2. Grid Refinement Study of Torczynski (1993) 

In Torczynski (1993), the Re = pu2h/n is based on upstream conditions. The variable u is 
the spatial average of the horizontal velocity u over h. The geometry is specified to match that 
of Kaiktsis et al. (1991). The upstream channel height h and step height a have values of h = 1 
and s = 0.94231, yielding a downstream channel height of H = 1.94231. The comer of the step 
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is at (*, y) = (1, 0). The channel extends a distance l—l upstream from the step and a distance 
L = 34 downstream from the step to preclude undue influence of the finite channel length on the 
flow at Re = 800. The following conditions are applied on the boundaries of the computational 
domain: u = v = 0 on the upper and lower channel walls, —p + pdu/dn — 0 and dv/ dn = 0 
on the outflow boundary, and u = [tanh(f/16)]uB(y) + [1 — tanh(</16)]up(y) and v = 0 on the 
inflow boundary and the step surface. Here, u#(y) = ma*[0, 3y(l — y)] is the correct boundary 
condition for flow over a backward-facing step and up(y) = 3(1 — y)(< 4- y)/(l + *)* is 
the Poiseuille flow observed infinitely far downstream whenever steady flow is asymptotically 
obtained. The initial velocity field is set equal to ti = tip(y) and t> = 0 throughout the domain. 
Here t> is the vertical velocity and p is the pressure. Thus, the above combination of boundary 
and initial conditions initially allows flow through the step surface so that the simulations can be 
initialized using an exact divergence-free solution of the Navier-Stokes equations. Furthermore, 
since the inflow boundary condition is varied smoothly in time from Poiseuille flow to flow 
over a backward-facing step, the flow experiences an order-unity transient that is probably 
strong enough to excite sustained unsteady behavior, if that is the appropriate asymptotic state 
for the numerical solution. 

The simulations were performed using the commercial code NEKTON v2.8, which employs 
a time-accurate spectral-element method with the Uzawa formulation (Nektonics, 1991). Let D 
be the dimensionality. Each element has N D velocity nodes located at Gauss-Lobatto Legendre 
collocation points, some of which are on the element boundaries, and {N — 2) D pressure nodes 
located at Gauss Legendre collocation points, all of which are internal. Within each element, 
the velocity components and the pressure are represented by sums of D-dimensional products 
of Lagrangian-interpolant polynomials based on nodal values. This representation results in 
continuous velocity components but discontinuous pressure at element boundaries. Henceforth, 
the quantity N is referred to as the element order, even though the order of the polynomials 
used to represent the velocity is N — 1. NEKTON employs mixed explicit and implicit temporal 
discretizations. To avoid solving a nonlinear nonsymmetric system of equations at each time 
step, the convective term is advanced explicitly in time using a third-order Adams-Bashforth 
scheme. All other terms are treated implicitly (backward Euler for the pressure and for the 
viscous terms). 

Three spectral-element grids of differing resolution, denoted L (low), M (medium), and H 
(high), are employed. Figure 5.8 shows the computational domain and the grid distribution 
of the three spectral element grids in which the distribution of nodes within each spectral 
element is not shown. The L grid with N = 9 is identical to the grid of Kaiktsis et al. 
(1991). Tables 5.1 and 5.2 show the results for a fixed time step of At = 0.10 at t = 800. 
This time step appears to satisfy the Courant restriction of the explicit portion of the temporal 
discretization. Four general classes of behavior are observed for the numerical solutions. First, 
“steady monotonic” denotes evolution of the numerical solution toward an asymptotically 
steady state. Second, “steady oscillatory” denotes evolution toward an asymptotically steady 
state with a decaying oscillation superimposed on the monotonic decay. Third, unsteady 
chaotic” denotes irregular transient behavior of the numerical solution that shows no indication 
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of evolving toward steady behavior. Fourth, “diverge” denotes a numerical solution terminated 
by a floating-point exception. In Tables 5.1 and 5.2, the first character under the “case” 
columns denotes the grid resolution £, M or H, the first number indicates the Reynolds number 
500 or 800 and the last two numbers indicate the order of the spectral element being used. For 
example, £807 means Re = 800 using the £ grid with N = 7. 

The extensive grid refinement study of Torczynski resulted in grid-independent steady-state 
numerical solutions for both Re = 500 and Re = 800. As the grid resolution is reduced below 
the level required to obtain grid independent solutions, chaotic-like temporal behavior occurred. 
The degree of grid resolution required to obtain a grid-independent solution was observed to 
increase as the Reynolds number is increased. 

5.4.3. Temporal Refinement Studies Using Knowledge from Dynamical Systems Theory 

The purpose of the present study is not to reaffirm Torczynski’ s study, or to target the spectral 
element method, or the work of Kaiktsis et al. ( 1 99 1 ), but rather to allude to how easily, without 
both temporal and spatial refinement study and knowledge of dynamical systems theory, similar 
erroneous conclusions could possibly arise in other CFD practices. As stated before, inaccuracy 
associated with underresolved grids (e.g., disappearance of the fine details of the flow due to 
large numerical diffusion and/or low order of the scheme) and linear instability associated with 
time discretizations are easier to detect than spurious numerical solutions. With the knowledge 
of possible nonlinear behavior of numerical schemes such as long time transients before a steady 
state is reached, numerically induced chaotic transients, numerically induced or suppressed 
chaos, existence of spurious steady states and asymptotes, and the intimate relationship among 
initial data, time step and grid spacing observed in discrete dynamical systems theory, we need 
to examine these cases in more detail. They can serve to illustrate the connection between the 
spurious numerics phenomena observed in simple nonlinear models and CFD computations. 

In the present study, in addition to grid refinement studies, temporal refinements are made 
on all of the underresolved grid cases in Tables 5.1 and 5.2 to determine if these cases sustain 
the same temporal behavior at a much later time or evolve into a different type of behavior. 
At t = 800, cases £506, £507, £508, £509, £811, M807 and M808 either exhibit “unsteady 
chaotic” or “steady oscillatory” behavior. We integrate these cases to t — 2000 to determine 
if a change in solution behavior occurs. From the phenomena observed in Keefe (1996) and 
others, t = 2000 might not be long enough for a long time transient or long chaotic transient 
to die out. There is also the potential of evolving into a different type of spurious or divergent 
behavior at a much later time. However, for this study it appears that t = 2000 is sufficient. For 
Re = 500, we also recomputed some of these cases with a sequence of At that bracketed the 
benchmark study of Torczynski. The A< values are 0.02,0.05,0.10,0.125,0.2,0.3,0.4, and 
0.5 for Re = 500. The CFL number for all of these cases is above 1 for At > 0.10. The reason 
for the investigation of At = 0.3, 0.4 and 0.5 is to find out, after the transients have died out, if 
the solution converges to the correct steady state for At that are a few times larger than 0.10. 

For Re = 800, we integrate £811 and Jl/808 with At = 0.10 and M807 with At = 
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0.02,0.05 and 0.10 to t = 2000. Aside from integrating to t = 2000, five different initial data 
were examined for cases il£807, AI809 and M 811 for A t = 0.10 to determine the influence of 
the initial data and the grid resolution on the final numerical solution. The five initial data are: 

(a) Uniform: u y v = 0 

(b) Shear layer: u = tz B (y) = ma*[0,3y(l — y)],v — 0 

(c) Solution from solving the steady Stokes equation (with no convection terms) 

(d) Torczynski (1993): u = i ip(y) = 3(1 - y)(< + y)/(l + = 0 

(e) Channel flow both upstream & downstream of step: Same as (d) except the boundary 

conditions 

The boundary conditions for (a), (b), (c) and (e) were parabolic inflow and no-slip at 
walls, whereas the boundary conditions for (d) were those of Torczynski (1993): u = 
[tanh(*/16)]u B (y) + [1 - tanh(t/16)]«p(y) and t» = 0. The CPU required to run the above 
cases ranged from less than a day to several days on a Sparc Center 2000 using one processor. 

The solution behaviors reported in Tables 5.1 and 5.2, with additional refinement study, now 
become Tables 5.3 and 5.4. Note that the “steady monotonic” cases in Tables 5.3 and 5.4 
are at t = 800 and the rest are at t = 2000 (if a divergent solution has not occurred earlier). 
The chaotic-like behavior evolves into a time-periodic solution beyond t = 800 for £506 and 
£507, whereas the chaotic-like behavior evolves into a time-periodic solution beyond t = 800 
for £811 and a divergent solution for M807. The “steady oscillatory” case £508 is slowly 
evolving to the correct steady state with an amplitude of oscillation of 10 -5 . The oscillation is 
not detectable within the plotting accuracy. The “steady oscillatory” time evolution of M 808 
is similar to that of £508. The numerical solutions with “steady oscillatory” and “steady 
monotonic” behavior at early stages of the time integration are almost identical at later stages 
of the time integration. They all converge to the correct steady state. The initial data study at 
Re = 800 with At = 0.10 is summarized in Table 5.5. It illustrates the intimate relationship 
between initial data and grid resolution. 

Figure 5.9 shows the streamlines for £509 (steady state solution) and £507 (spurious 
time-periodic solution). Figure 5.10 shows the streamlines for J5T809 (steady solution) and £811 
(spurious time-periodic solution) and the corresponding grids with the distribution of the nodes 
of the spectral elements shown. Note that even for the £ grid using N — 11 (£811), the grid 
spacings are very fine and yet a spurious time-periodic solution was obtained. 

Figures 5.1 1 - 5.14 show the vertical velocity time histories at (as, y) = (30, 0) advanced to 
a time of t = 2000 of selected runs for both Reynolds numbers and various At. They illustrate 
the different spurious behaviors of the underresolved grid cases at t < 2000. The amplitude of 
the spurious time-periodic solutions remains uniform for £506 for 0.02 < At < 0.2 (Fig. 5.1 1) 
but not for £507 (Fig. 5.12). A counter-intuitive behavior was observed for the £507 case. For 
£507, the amplitude of the spurious time-periodic solution remains constant for At = 0.02 and 
0.05 but decreases to a significantly lower value for the large time-step range. One would expect 
the opposite effect on the height of the amplitude. In addition, the two distinct amplitudes of the 
periodic solution indicate the existence of two finite ranges of At where the numerical solutions 
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converge to two distinct spurious time-periodic solutions for £507. 

Figure 5.14 shows the vertical velocity time histories at (*,y) = (30,0) for M 807 with 
At = 0.02, 0.05 and 0.10, and 1811 for At = 0.10. Case M807 diverges at t = 1909.2 for 
At = 0.02, at t = 972.4 for At = 0.05, and at t = 827.77 for At = 0.10. The time histories 
for these three time steps appear to show chaotic-like behavior if one stops the computations at 
t = 800. The bottom plot of Fig. 5.14 shows the vertical velocity time histories advanced to a 
time of t = 2000 for £811 with At = 0.10. It shows the definite time-periodic spurious solution 
pattern. On the other hand, the time history for this case appears to show an aperiodic-like 
pattern if one stops the computation at t = 800. Note that the £809 grid case was used by 
Kaiktsis et al. (1991) and they concluded that transition has occurred at Re = 800. 

In summary, without the temporal refinement study (longer time integration), the £506, 
£507, £811 and Af807 cases can be mistaken to be chaotic-like (or aperiodic-like) flow. 
Although the time history up to t = 800 appears chaotic-like, one cannot conclude it is chaotic 
without longer transient computations. Comparing Tables 5.1 and 5.2 with 5.3, 5.4 and 5.5, 
one can conclude that with transient computations that arc 2.5 times longer than Torczynski’s 
original computations, what appeared to be aperiodic-like or chaotic-like behavior at earlier 
times evolved toward either a time-periodic or divergent solution at later times. These temporal 
behaviors appear to be long time aperiodic-like transients or numerically induced chaotic-like 
transients. For Re = 800, five different initial data were examined to determine if the 
flow exhibits strong dependence on initial data and grid resolution. Results showed that the 
numerical solutions are sensitive to these five initial data. The refinement study also revealed a 
nonstandard guideline in grid clustering or grid adaptation. Traditional grid refinement and grid 
adaptation methods concentrate on regions with strong gradients, shock waves, slip surfaces 
and fine structure of the flow, and de-refine regions of smooth flows. As can be seen in Figs. 
5.9 and 5.10, the flow down stream of the backward facing step is very smooth, yet a fine 
grid is needed in this region in order to obtain the correct numerical solution. It is postulated 
that proper nonlinear wave propagation is hampered by the underresolved grid. This behavior 
may be related to spurious discrete traveling waves. A separate investigation is needed and 
is beyond the scope of the present paper. See Yee et al. (1991), Griffiths et al. (1992b) for 
a discussion. A suggestion to minimize spurious asymptotes is discussed in Yee & Sweby 
(1995b, 1996b). Note that the results presented pertain to the characteristic of the studied 
scheme and the DNS computations. However, if one is certain that Re = 800 is a stable 
steady flow, a non-time-accurate method such as time-marching to obtaining the steady-state 
numerical solution would be a more efficient numerical procedure. 

5.5. Spurious Behavior of Time-Lag Coupling of a Fluid-Structure Interaction 

This section discusses the spurious behavior of a semi-implicit method that is commonly used 
in combustion and fluid-structure interaction analysis. These types of problems are commonly 
mathematically stiff and highly nonlinear, and consist of a large number of strongly-coupled 
equations. The simulation considered employs the unsteady 2-D compressible Navier Stokes 
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equations coupled with a two-equation elastic model on overlapping grids involving stationary 
and deforming meshes. 

A common practice in time-accurate aeroelastic computations is to utilize well validated 
implicit Navier-Stokes algorithms that were developed for complex flowfields over 3-D 
nondeforming bodies and extend them to include aeroelastic effects. The simplest approach 
to extend these algorithms is to lag the effects of moving/deforming structures by one time 
step (Smith 1989, Guruswamy 1990, Morton & Beran 1995), allowing current algorithms to 
be used in updating the aerodynamic variables. After the aerodynamic loads are determined, a 
structural module is called to update the position and shape of the body. A disadvantage of this 
semi-implicit strategy is the fact that regardless of the temporal accuracy of the aerodynamic 
and structural algorithms, the loose coupling introduces an 0( At) error, and may impose more 
stringent stability criteria. . Besides inlroducing^undesirable. phase lag, it may also introduce 
spurious solutions as a result of the semi-implicit procedure. 

5.5.1. Elastically Mounted Cylinder 

An elastically mounted cylinder model problem is used to demonstrate the occurrence of 
spurious solutions when using the lagged structures approach. A circular cylinder is mounted 
in a freestream of velocity V*, with linear springs in both coordinate directions as depicted in 
Fig. 5.15. The aeroelastic cylinder is an attractive model problem for two reasons. First, it 
displays nonlinear unsteady flowfield physics associated with separation and vortex shedding, 
and, secondly, there are both numerical and experimental data available for comparison (Alonso 
et al. 1993, Blackburn & Kamiadakis 1993). 

The governing equations used to model the aerodynamic system are the compressible 
laminar Navier-Stokes equations. The governing equations for the two-degree-of-freedom 
model (Alonso et al. 1995, Blackburn & Kamiadakis 1993) in dimensional form are 

mx ta + Ci'a + Kx ta = D, (5.18) 

my ea + Cy ta + Ky ta = X, (5.19) 

where m, C, K , D, and X are the mass, coefficient of structural damping, coefficient of spring 
stiffness, drag, and lift per unit span, respectively, and x ea and y ea are the horizontal and 
vertical positions of the center of the cylinder. 

The governing equations are solved on a deforming mesh overlapping a stationary mesh with 
a Beam-Warming approximate factored algorithm modified to include Newton-like subiteration 
of index p, coupled with an ordinary differential equation structural solver, also in subiteration 
form. The temporal discretization is either the backward Euler (first-order) or the three-level 
backward differentiation (3-level BDF, second-order) and is linearized about the solution at 
subiteration level p. The spatial derivatives of the Navier-Stokes equations are approximated 
by second-order central differences and common forms of both implicit and explicit nonlinear 
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dissipation. With a sufficient number of subiterations, this approach becomes a fully-implicit 
first- or second-order accurate aeroelasticity solver. All solutions of this work were computed 
using the backward Euler temporal discretization without subiteration to model a lagged 
structures approach. The solutions were then compared with the solutions of the fully-implicit 
approach. Details of the solver can be found in Morton et al. (1997). Morton et al. illustrate 
the importance of using the fully-coupled fully-implicit second-order approach for the fluid- 
structure interaction. Surface boundary conditions are comprised of no slip, adiabatic wall, and 
the inviscid normal momentum equation. Freestream conditions are specified along the outer 
boundary inflow and extrapolation in the horizontal coordinate is implemented at the outer 
boundary outflow. Periodic conditions are applied along the overlap boundary (due to the O 
grid topology). The method’s accuracy was verified through comparison with numerical and 
experimental solutions reported in Morton et al. (1997). 

5.5.2. Numerical Results 

Elastically mounted cylinder solutions were computed for a variety of time steps with the 
lagged structures/no subiteration approach. The baseline structural parameters used for all cases 
were 


Re = 500, Moo = 0.2, ( = 1, fi. = 5, tt = 4 (5.20) 

where ( is the nondimensional structural damping coefficient, ft, is the mass ratio and u is 
the reduced velocity. See Morton et al. (1997) for details of the physical parameters. The 
computational grid has 384 evenly spaced points around the cylinder, and 96 points in the radial 
direction. A nondimensional spacing of 0.0005 was specified normal to the surface, and the 
grid was geometrically stretched to a maximum radius of 50 cylinder diameters. 

The fluid-structure interaction system was initialized by computing a static cylinder time- 
periodic solution with a time step of At = 0.01. Once this solution was determined to be 
periodic, the cylinder was allowed to move in both coordinate directions in response to the 
periodic shedding of vortices. The cylinder established a new periodic solution characterized 
by oscillations in the * and y coordinates of the cylinder center. This initial periodic solution 
was then used to compute a set of solutions for increasing and decreasing time steps. 

Refinement in time step produced an asymptotic solution with a nondimensional frequency 
(Strouhal number) of St = 0.2256. Solutions for the most refined time step up to a time step 
of 0.02 were sinusoidal with a single frequency and amplitude. As the time step was increased 
in this range, the amplitude of vertical motion increased and the frequency of oscillation 
decreased monotonically. Solutions for time steps greater than 0.02 have additional frequency 
content not evident in the smaller time step solutions. To determine the frequency content, a 
power-spectral-density (PSD) analysis was performed with MATLAB. 

Figures 5.16 and 5.17 show the cylinder center vertical motion time histories and the 
corresponding PSD analysis for six different time steps (0.01 < At < 0.06 with an increment 
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of .01). The length of the total time integrations, determined by monitoring the time evolution 
solution behavior, increases as Af increases. The total time integrations for the six time steps are 
t = 380, 900, 1400, 1800, 2200 and 2600, respectively. The lengths of these time integrations 
are also guided by the knowledge gained from previous sections and Yee & Sweby (1996a,b). It 
is interesting to see the various spurious behavior as a function of At. The form of the solutions 
evolves from a sinusoidal periodic solution to periodic solutions with more than one frequency, 
and eventually to aperiodic chaotic-like patterns. 

Figures 5.17a,b show the single frequency associated with time steps of 0.01 and 0.02 with 
a trend toward lower frequency with increasing time step. The amplitudes of the sinusoidal 
motion are affected by the size of the At. The PSD analysis for the At = 0.03 solution (Fig. 
5.17c) shows the same trend for the dominant frequency but an additional lower frequency 
component is evident. This .additional spurious frequency, is responsible for the aperiodic- like 
motion of the cylinder (Fig. 5.16c). 

It is interesting to note the change in character of the solution for At = 0.04 (Fig. 5.16d). 
The solution is time-periodic with several distinct local minima and maxima within one cycle. 
The PSD analysis shows three distinct frequencies (Fig. 5.17d) with the two spurious low 
frequencies dominating the frequency closest to the asymptotic solution frequency. 

The solutions for At = 0.05 and 0.06 (Fig. 5.16e,f) are chaotic-like. The PSD analyses 
(Fig. 5.17e,f) for both time steps show a spectrum of frequencies as opposed to the few distinct 
frequencies seen in lower At solutions. The dominant frequencies are similar for both time step 
solutions with the At = 0.06 solution showing very large PSD in the low frequency range. 

These time histories indicate a counter-intuitive behavior. The solution changes from a 
periodic pattern for At = 0.02 to an aperiodic pattern for At = 0.03, and then back to a 
periodic pattern for A* = 0.04 before the onset of chaotic-like behavior for At = 0.05 and 
0.06. 

Solutions were computed with the fully-coupled approach with subiterations to determine if 
spurious solutions were evident at the same time steps. Figure 5.18 depicts a comparison of 
fully-implicit versus lagged structures solutions for At = 0.01 and 0.05. The spurious behavior 
was not evident for At = 0.05 for the fully-implicit approach. In addition to the spurious 
behavior manifested by the lagged structures for slightly larger time steps, smaller time steps, 
although producing the correct solution, exhibit a time lag over the fully-coupled case. 

In summary, for time steps greater than 0.02, the model exhibits spurious solutions when 
the loosely-coupled implicit approach is employed. In some cases the numerical solutions were 
not chaotic but were still spurious and time-periodic, making it difficult for the researcher to 
determine if the solution is representative of the true physics of the problem. Fortunately, a 
fully-implicit structural coupling eliminated the spurious solutions for time steps much greater 
than those associated with the spurious lagged-structures solutions for the given model problem. 
Large scale full aircraft computations are expensive and therefore it is tempting for researchers 
to trade efficiency for accuracy by increasing the time step. This may lead to spurious solutions 
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that are difficult to detect without a comprehensive time step refinement study. 


5.6. Strong Dependence on Initial Data & Underresolved Grids of a 3-D Simulations of 
Vortex Breakdown on Delta Wings 

This section discusses the spurious behavior of underresolved grids observed in Visbal 
(1995a,b, 1996) using the backward Euler fully-implicit temporal discretization for a 3-D 
vortex breakdown on delta wings. For certain initial angles of attack, spurious time-periodic 
and spurious chaotic-like temporal behavior occurs as the grid resolution is reduced. The 
coarse grid used is actually finer than grids commonly used in full aircraft simulations. In view 
of the fact that some experimental studies have reported the existence of vortex breakdown 
static hysteresis on delta wings and others have not, the coarse grid case could be mistaken to 
exhibit similar nonunique solutions behavior to that of some of the experimental data if a grid 
refinement study is not made. 

5.6.1. Flow Configuration, Governing Equations and Numerical Procedure 

The vortical flows encountered by agile aircraft at high-angle-of-attack exhibit a variety of 
complex, nonlinear aerodynamic phenomena not yet fully understood. Among these is ‘vortex 
breakdown’ or ‘vortex bursting’ which represents a sudden disruption of the well-organized 
leading-edge vortex present above slender wings at high incidence. Vortex breakdown is 
typically characterized by reverse axial flow in the vortex core, and by marked flow fluctuations 
downstream of the breakdown inception location. The sudden onset of vortex breakdown 
and its effects on aerodynamic loads severely impact aircraft stability and control and may 
result in reduction of the operational envelope of high-performance aircraft. In addition, 
the breakdown induced fluctuations may promote undesirable fluid/structure interactions on 
aircraft components intersecting the vortex path. Such is the case of tail-buffet present in twin- 
tailed aircraft. Reviews of important aspects of vortex breakdown on delta wings have been 
provided in Lee & Ho (1990), Rockwell (1993) and Visbal (1995a). These indicate that despite 
recent progress, vortex breakdown still remains a challenge in its fundamental understanding, 
prediction and control. Further insight into this flow phenomenon could be achieved by 
systematic experimental and computational studies describing the complex three-dimensional, 
unsteady structure of vortex breakdown flow fields. The present section describes spurious 
solutions encountered while performing a detailed computational study in Visbal (1995a,b, 
1996) of the spiral vortex breakdown structure above a slender delta wing. 

The flow configuration considered in Visbal (1995a,b, 1996) consists of a flat-plate delta 
wing with a sweep angle of A = 75° and a freestream velocity £/«, depicted in Fig. 5.19. 
The freestream Mach number and the Reynolds number based on the centerline chord are 0.2 
and 9.2 x 10*, respectively. The sweep angle and Reynolds number were selected to permit 
comparison with the extensive experiments of Magness (1991). 

The governing equations for the present simulations are the unsteady, 3-D compressible 



83 


full Navier-Stokes equations supplemented by the perfect gas law, Sutherland’s viscosity 
formula and the assumption of a constant Prandtl number (Pr = 0.72). These equations are 
written in strong conservation-law form using a general coordinate transformation. For the low 
Reynolds number considered, the flows are assumed to be laminar, and no turbulence models 
are employed. 

The governing equations are numerically solved employing the implicit, approximate- 
factorization, Beam-Warming algorithm (Beam & Warming 1978). The scheme is formulated 
using backward Euler time-differencing and second-order central difference approximations for 
all spatial derivatives. Fourth-order non-linear dissipation terms are added to control odd-even 
decoupling. Newton-like subiterations are also incorporated in order to reduce linearization and 
factorization errors, thereby improving the temporal accuracy and stability properties of the 
algorithm. A vectorized, time-accurate, three- dime nsional-Navier- Stokes solver (FDL3D1) has 
been developed to implement the previous scheme. This code has been validated extensively 
for both steady and unsteady flow fields (see Visbal 1996 and references therein). 

5.6.2. Grid Structure and Boundary Conditions 

The computational grid topology for the flat-plate delta wing is of the H-H type and is 
obtained using simple algebraic techniques. For this mesh topology, the boundary conditions 
are implemented in the following manner. On the lower, upper, lateral and upstream boundaries, 
characteristic boundary conditions are specified. On the downstream boundary, through which 
the vortex exits, flow variables are extrapolated from the interior. It should be noted that 
since the grid is smoothly stretched toward the downstream boundary, the vortical structure has 
almost entirely dissipated prior to reaching the end of the computational domain. On the wing 
surface, no-slip, adiabatic conditions are applied in conjunction with the usual zero normal 
pressure gradient approximation. In this study, only half of the delta wing is considered and 
symmetry conditions are imposed along the mid-plane of the wing. This is done in order to 
provide better numerical resolution of the spiral breakdown with a given number of grid points, 
at the expense of not being able to resolve asymmetric effects. This approach is therefore 
limited to angles of attack for which breakdown location is not too far upstream as to result in 
changes in breakdown structure due to the interaction of the two vortices. 

In order to assess numerical resolution effects, three different grid sizes, denoted as Grid 1 , 
2 and 3 have been employed with streamwise (X), spanwise (F), and normal directions (Z) 
respectively (see Fig. 5.19): 

Gridl: 98 x 115 x 102 

Grid 2: 159 x 107 x 149 

Grid 3: 209 x 107 x 149 

For Grid 1, the streamwise spacing over the wing is AX / C = 0.02, where C denotes the wing 
centerline chord. The spanwise and normal spacing near the vortex axis at the trailing edge 
(X/C = 1.0) are A Y/C = 0.007 and AZ/C = 0.008 respectively. In Grid 2, the streamwise 
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spacing is halved (AX/C = 0.01) in the region where vortex breakdown is anticipated 
(0.5 < X/C < 1.2). In addition, both the spanwise and normal grid spacings in the vortex axis 
are reduced by a factor of two. Finally, the finest mesh Grid 3 is obtained from Grid 2 by further 
decreasing the streamwise spacing in the breakdown region to A X/C = 0.005, bringing the 
cell aspect ratio in the vortex core closer to one. From a grid resolution study presented in 
Visbal (1996), it was concluded that Grid 2 was sufficient to capture the basic structure of a 
spiral vortex breakdown. The results to be described in the next section pertain to Grids 1 and 
2 only. 


5.6.3. Results 

According to the experiments of Magness (1991), vortex breakdown moves upstream over 
the delta wing when the angle of attack is slowly increased (i.e. in a quasistatic manner) beyond 
a critical value a cr ~ 30.7°. In order to study computationally this quasistatic onset of vortex 
breakdown, calculations were performed initially near a er using a 1° increment in angle of 
attack and employing two of the grids previously described. Although smaller steps in incidence 
would have been desirable, this was not computationally feasible. The solutions obtained prior 
to breakdown onset (a = 30°) using Grid 1 (Fig. 5.20a) and Grid 2 (Fig. 5.20d) are found to be 
in reasonable agreement with each other. The computed onset of vortex breakdown occurred for 
Grid 1 when a was increased from 31° to 32°, whereas for Grid 2, it took place between 30° and 
31°, in closer agreement with experiment (Magness, 1991). The process by which breakdown 
appeared in the near-wake and moved over the wing was also found to be qualitatively similar 
for both grids. Based on these comparisons, one would conclude that the effect of numerical 
resolution on the quasistatic breakdown onset is small. However, as described below, the 
non-linear dynamic behavior near a cr was found to be affected significantly by numerical grid 
resolution when large increments in the initial angle of attack are imposed. 

Initial Angles of Attack: As discussed in Visbal (1995a), some experimental studies have 
reported the existence of vortex breakdown static hysteresis on delta wings (i.e. multiple, time- 
asymptotic solutions for a given static angle of attack). Motivated by these findings, the existence 
of non-unique solutions for different initial conditions was investigated numerically. However, 
instead of observing static hysteresis phenomena similar to that in some of the experimental 
studies, spurious behavior due to the numerics was encountered. Computations were performed 
on Grids 1 and 2 for a = 30° using three different jump-start initial angles of attack a 0 . 
With the exception of the initial conditions and grids, all remaining numerical parameters 
(i.e., time step, damping coefficient, time-marching procedure and boundary conditions) were 
kept the same. In all of the computations, the fixed non-dimensional At + = 0.0005. This 
non-dimensional Af + is equal to the dimensional At times Uoo/C. The choice of the fixed At 
is based on the study reported in Visbal (1996). The three initial angles of attack ao denoted by 
Cases a,b and c are: 

Case a: ao = 29° 

Case b: ao = 25° 
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Case c : a 0 = 17° 

Grid 1: The computed lift coefficient histories obtained on Grid 1 for the above three initial 
conditions are shown in Fig. 5.21 which exhibit three distinct numerical solutions. The Case a 
solution corresponds to the columnar (no breakdown) solution shown in Fig. 5.20a. For Case 
b with a mild jump in the initial angle of attack ao = 25° to the desired a = 30°, the lift 
coefficient history appears to be time- periodic. This solution asymptotes to a flow containing a 
mild vortex breakdown near the wing trailing edge as shown in Fig. 5.20b. The frequency spectra 
of the streamwise velocity fluctuations at a point within the reverse flow region in the vortex 
core is given in Fig. 5.22. It indicates the nearly periodic character of this spiral breakdown. For 
Case c, a third distinct solution was obtained for a large jump in the initial angle of attack from 
a 0 = 17° to the desired a = 30°. The lift coefficient time history appears to be chaotic-like. 
This flow field exhibits breakdown well upstream of the trailing edge, as shown in Fig. 5.20c. 
Also, the corresponding velocity fluctuations associated with this stronger breakdown display 
multiple frequencies (Fig. 5.23). Both Cases b and c were run till = tUoo/C ~ 80, during 
which time no tendency was observed for the breakdown to leave the wing. For the purpose of 
comparison with breakdown computations in tubes or isolated vortices, < + = 80.0 corresponds 
in the present case to 2400 characteristic times based on the vortex core radius at Xj C = 0.5. 
Comparison of the three solutions in terms of the X— component of vorticity in the vortex core 
at a location upstream of breakdown ( X/C = 0.4) is shown in Fig. 5.24. All flow solutions are 
seen to be essentially the same at this upstream location, indicating that the different cases are 
not simply due to lags in the development of the vortex following a change in angle of attack. 

Grid 2: In order to investigate the validity of the qualitative behavior displayed by the 
finite-difference solutions computed on Grid 1, calculations were also performed on Grid 2 for 
a = 30° using the same three initial conditions. The corresponding lift coefficient histories are 
shown in Fig. 5.25. Instead of obtaining three distinct solutions as in the Grid 1 case, the lift 
coefficient time histories for the three initial data using Grid 2 evolved to a similar behavior at 
a later stage of the time integration. The flow structure of Case a is given in Fig. 5.20d. On 
this finer mesh, the solution exhibits unsteady boundary-layer separation near the wing trailing 
edge which results in the observed lift coefficient fluctuations. Notice, however, that the mean 
Cl is in close agreement with the corresponding value on Grid 1 (Fig. 5.21, a 0 = 29°). This 
unsteady phenomenon, captured with the improved spatial resolution, is entirely separate from 
vortex breakdown and does not significantly affect the vortex core (as seen by comparison of 
Figs. 5.20a and 5.20d). For Case b, vortex breakdown appears transiently above the wing and 
penetrates upstream up to X/C ~ 0.85, as described in Visbal (1995b). However, later on 
breakdown moves downstream off the wing and into the wake, albeit at a very slow rate. By 
t+ » 30.0, the computed lift coefficient for Case b has reached the same level and characteristic 
behavior found for Case a, and the flow structure (not shown) is essentially the same as that 
of Fig. 5.20d. For Case c, a strong vortex breakdown is induced upstream of X/C = 0.8 (see 
Visbal 1995b). As before, it proceeds downstream at a very low speed, eventually leaving the 
wing. The computed lift coefficient (Fig. 5.25) is again observed to increase toward the same 
level attained in Cases a and b. Examination of the flow fields above the wing for Cases a, b 
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and c in Grid 2 indicated that they all asymptote to the same columnar solution. 

The previous results clearly indicate the different dynamic qualitative behavior displayed by 
the discrete equations on the two grids employed. With the coarser mesh, non-unique solutions 
which exhibit vortex breakdown are achieved for a < a cr if a sufficiently large A a is imposed 
in the initial conditions. When the finer grid is utilized, however, the non-unique solutions are 
no longer found. Vortex breakdown is still induced transiently above the wing but it eventually 
moves into the wake, albeit at very low speed. The behavior computed on Grid 2 is in agreement 
with the experiments of Magness (1991). In these experiments it is found that when the wing is 
pitched at a high rate from a = 5° to a = 30° (below a cr ), breakdown occurs over the wing 
up to X/ C « 0.6. Upon cessation of the motion, however, breakdown leaves the wing, and no 
static hysteresis (due to multiple solutions) was found. 

5.6.4. Implications of Spurious Behavior 

The above discrepancies in qualitative behavior obtained with the two levels of spatial 
resolution indicate either that (1) spurious time-asymptotic solutions containing breakdown can 
be obtained on underresolved grids, or that (2) the basins of attractions (allowable initial data) 
of the various solutions (if they exist) change vastly with mesh spacing as observed in simple 
model problems (Yee & Sweby 1996) and in Sections II and III. Regardless of the cause, we 
can safely conclude that for the given model problem and flow field conditions, the nonunique 
solutions exhibited by Grid 1 are numerical artifacts. Although a precise explanation cannot be 
offered at present, the first possibility is postulated as being the reason for these discrepancies. 
It is possible that the larger mesh spacing in Grid 1 provides an artificial mechanism for 
wave-trapping if breakdown is transiently induced. However, further work is clearly required 
to verify this hypothesis. In the meantime, these results serve to point out the danger associated 
with interpreting complex flow behavior computed with underresolved meshes. Since Grid 1 
is actually finer than grids commonly used in full aircraft simulations, the present results have 
important implications for practical CFD studies where systematic resolution assessments are 
not always feasible. The spurious dynamic behavior discussed is particularly relevant, since 
high-angle-of-attack calculations are typically started (in order to minimize computer time) 
with abrupt uniform-flow initial conditions or with discrete jumps in angle of attack which may 
induce transient vortex breakdown. 

VI. Concluding Remarks 

The need for the study of dynamics of numerics is prompted by the fact that the type 
of problem studied using CFD has changed dramatically over the past decade. CFD is also 
undergoing an important transition, and it is increasingly used in nontraditional areas. But even 
within its field, many algorithms widely used in practical CFD applications were originally 
designed for much simpler problems, such as perfect or ideal gas flows. As can be seen in 
the literature, a straightforward application of these numerical methods to high speed flows, 
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nonequilibrium flows, advanced turbulence modeling or combustion related problems can lead 
to wrong results, slow convergence, or even nonconvergent solutions. The need for new 
algorithms and/or modification and improvement to existing numerical methods in order to 
deal with emerging disciplines is evident. We believe the nonlinear dynamical approach for 
CFD can contribute to the success of this goal. The first step toward achieving this goal is 
to understand the nonlinear behavior, limits and barriers, and to isolate spurious behavior of 
existing numerical schemes. 

We have revealed some of the causes of spurious phenomena due to the numerics in an 
attempt to improve the understanding of the effects of numerical uncertainties in CFD. We have 
shown that guidelines developed using linearization methods are not always valid for nonlinear 
problems. We have gained an improved understanding of long time behavior of nonlinear 
problems and nonlinear stability, convergence and reliability of time-marching approaches. We 
have learned that numerics can introduce and suppress chaos and can also introduce chaotic 
transients. The danger of relying on numerical tests (e.g., direct numerical simulation) for the 
onset of turbulence and chaos is evident. 

We illustrated with practical CFD examples that exhibit properties and qualitative behavior 
similar to those of elementary examples in which the full dynamical behavior of the numerics 
can be analyzed. Some of these CFD examples were chosen based on their non-apparent 
spurious behaviors that were difficult to detect without extensive grid and temporal refinement 
studies and without some knowledge from dynamical systems theory. In all of the underresolved 
grid cases, the grids actually are finer than grids commonly used in realistic complex CFD 
simulations. The three semi-implicit methods considered in Section V are typical numerical 
procedures employed in active research areas such as chemically reacting flows, combustion, 
fluid-structure interactions, DNS and large eddy simulations (LES). These studies serve to point 
out the various possible dangers of misinterpreting numerical simulations of realistic complex 
flows that are constrained by available computing power. 

The observed spurious behavior related to underresolved grid cases is particularly relevant to 
DNS and LES. Spatial resolutions in DNS and LES are largely dictated by the computer power, 
especially when numerical algorithms other than high accuracy spectral methods are employed. 
LES, by design, filters out the small scales from the nonlinear Navier-Stokes equations. The 
effect of small scales on the large scale motion is accounted for by the subgrid scale model. 
Spurious behaviors due to underresolved grids in LES can play a major but different role from 
DNS. A dynamical approach study on the effect of underresolved grids is postulated to be 
useful in pinpointing the limitation of DNS and LES approaches and their associated spurious 
behavior that might be otherwise difficult to detect. These are subjects of future investigations. 

As can be seen, recent advances in dynamics of numerics showed the advantage of adaptive 
step size error control for long time integration of nonlinear ODEs. Although much research 
is needed to construct suitable yet practical similar adaptive methods for PDEs, these early 
developments lead our way to future theories. There remains the challenge of constructing 
adaptive step size control methods that are suitable yet practical for time marching to the 
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steady states for aeronautical CFD applications. Another even more challenging area is the 
quest for an adaptive numerical scheme that leads to guaranteed and rapid convergence to the 
correct steady-state numerical solutions. These two key challenges are particularly important 
for CFD. We conclude the paper with the following guidelines to minimize spurious dynamics 
in time-marching to the steady state. 

Some Guidelines to Minimize Spurious Dynamics'. Due to the spurious dynamics intro- 
duced by the numerics, one usually will not be able to map out the complete numerical basins 
of attraction and bifurcation diagrams for the entire problem in practical situations. Only in 
isolated situation with a particular physical problem and numerical method combination such 
as the one studied in Lafon & Yee (1991, 1992) are continuation methods able to locate all of 
the essential spurious branches of the bifurcation curves. 

On the other hand, continuation methods are widely used in dynamical systems when one 
wants to understand certain properties of key branches of the bifurcation curve, especially if 
one knows (or can ascertain by other means) a starting solution on that particular branch. For 
example, in the Taylor-Couette flow problem, extensive use is made of continuation methods 
to map out the critical Reynolds number when the flow behavior undergoes drastic changes in 
flow patterns, since in this case we know how the flow behaves for the low Reynolds number 
case. Another example is in Bailey and Beam (1991) where continuation methods were used 
to study the hysteresis behavior of the flow of an airfoil in terms of angles of attack for the 
steady PDEs. In this case, the flow behavior is readily obtained for low angles of attack 
(before hysteresis). Most of the use of the continuation method so far is focused on elliptic 
PDEs or steady PDEs. These studies seldom address the possibilities of spurious dynamics 
due to the numerics, especially for IBVPs using time-marching approaches. We note that a 
shortcoming in association with solving the steady PDEs is that a small radius of convergence or 
nonconvergence of the numerical solution is often encountered even with the aid of multigrid, 
preconditioners, and/or relaxation methods, especially when the PDEs are of the mixed type 
(e.g., the steady inviscid supersonic flow over a blunt body). In addition, the solutions obtained 
do not distinguish between stable and unstable steady states. 

Besides the study in Lafon & Yee (1991, 1992), here we propose a further step of applying 
this technique to the discretized counterparts of the time-dependent PDEs in order to avoid 
spurious asymptotes due to unknown initial data. The idea relies on the knowledge of a known 
or a reliable numerical solution on the correct (non-spurious) branch of the bifurcation curve 
as a function of the physical parameter of interest. The logic is that if one starts on the correct 
branch, one avoids getting trapped on any of the spurious branches. Also the issue of unknown 
initial data related to time-marching approaches is avoided or can be minimized. Details of the 
approach and numerical examples will be reported in a forthcoming paper. Here we will give a 
short narrative summary of the procedure. 

In many fluid problems the solution behavior is well known for certain values of the physical 
parameters but unknown for other values. For these other values of the parameters, the problem 
might become very stiff and/or strongly nonlinear, making the available numerical schemes (or 
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the scheme in use) intractable. In this situation, continuation methods in bifurcation theory can 
become very useful. If possible, one should start with the physical parameter of a known or 
reliable steady state (e.g., flow behavior is usually known for low angles of attack but not for 
high angles of attack). One can then use a continuation method such as the pseudo arclength 
continuation method of Keller (1977) (or the recent developments in this area) to solve for the 
bifurcation curve as a function of the physical parameter. The equations used are the discretized 
counterpart of the steady PDEs or the time-dependent PDEs. If time-marching approaches 
are used, a reliable steady-state numerical solution (as a starting value on the correct branch 
of the bifurcation curve for a particular value of the physical parameter) is assumed. This 
starting steady-state numerical solution is assumed to have the proper time step and initial 
data combination and to have the grid spacing fine enough to resolve the flow feature. The 
continuation method will produce a continuous spectrum of the numerical solutions as the 
underlying physical parameter is varied until it arrives at a critical value p c such that it either 
experiences a bifurcation point or fails to converge. Since we started on the correct branch of 
the bifurcation curve, the solution obtained before that p c should be more reliable than if one 
starts with the physical parameter in question and tries to stretch the limitation of the scheme. 
Note that by starting a reliable solution on the correct branch of the bifurcation curve, the 
dependence of the numerical solution on the initial data associated with time-marching methods 
can be avoided. 

Finally, when one is not sure of the numerical solution, the continuation method can be used 
to double check it. This approach can even reveal the true limitations of the existing scheme. 
In other words, the approach can reveal the critical physical parameter for which the numerical 
method breaks down. On the other hand, if one wants to find out the largest possible time step 
that one can use for a particular problem and physical parameter, one can also use continuation 
methods to trace out the bifurcation curve as a function of the time step. In this case, one can 
start with a small time step with the correct steady state and observe the critical time step as it 
undergoes instability or bifurcation. 
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Scheme 

Fixed points 

Stable range 

Explicit Euler 

i 

0 <r <2 

Modified Euler 

i 

0 <r <2 


l + 2 /r 

0 <r<-l + V5 = 1-236 


2/r 

2 < r < 1 + V5 = 3*236 

Improved Euler 

i 

0 < r <2 


[2 + r ± Vr 2 - 4] 

2 < r < V8 == 2-828 


2 r 

Heun 

1 

0< r < 1 + (Vl7 + 4) ,/3 - (Vl7 + 4)- ,/3 = 2-513 


* 

4-9137 <r <4-9552 


♦ 

6-4799 <r< 6-4853 


* 

6-74405 < r < 6-74575 

R-K 4 

1 

0<r<5 + (W + 5\/29)'' 3 



+ (^f - 5V^9 ) ,,j = 2-785 


* 

2-785 < r < 3-4156 


* 

2-746 < r < 3-456 


Table 3.1. Fixed points of Runge-Kutta methods for du/dt = au(l - u). 


Scheme 


Fixed points 

Stable range 

Explicit Euler 


1 

2 

0<r <8 

Modified Euler 


1 

2 

1(1 + Vl - 32/r] 
^[1 ± Vl - 8/r] 
i[3 - Vl - 32/r] 

0<r <8 

32 <r< 32-014067 
8< r <4(1 + V3) = 10-928 
32 <r <32 014067 

Improved Euler 


1 

1 ± Vl - 8/r 
2 

0<r<8 
8 < r < 12 


1 

2 

Vl - 12/r Vl+4/r 
4 4 

12 < r < 4(1 + V6) = 13-798 


1 

2 

Vl - 12/r Vl + 4/r 
+ 4 ± 4 

12 < r < 4(1 + V6) = 13-798 

Heun 


1 

2 

0 < r < 4(1 + (Vl7 + 4) 1/3 

-(VT7 + 4)- ,/3 )* 10-051 

R-KA 


\ 

2 

0 < r < t + 4(tt + 3 V29) 1 ' 3 

+ 4(^-5V29) l/3 = 11-14 


Table 3.2. Fixed points of Runge-Kutta methods for du/dt = au(l - u)(0.5 - u). 
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S(u) 

(3.11) 

(3.17) 

(3.18) 

(3.19) 

(3.21) 

“(1 -«) 

0 

2 

2 

6 

14 

u(i -«)(!-«) 

0 

6 

6 

24 

78 


Table 3.3. The number of spurious fixed points of Runge-Kutta methods 
for (3.7a) and (3.23). 


Equation 

Period 2 orbits 

Stable range 

u' = au( 1 — u) 

r + 2 ± Vr 2 — 4 

2r 

2<r<V6“ 2-4495 

u' = au(\ - u)Ci - u) 

1 ± Vl - 8/r 
2 

8 < r < 12 


1 Vl-12/r Vl+4/r 

2 4 ± 4 

12 < r < 14 


1 Vl-12/r Vl+4/r 

2 + 4 1 4 

12 < r < 14 


Table 3.4. Period 2 fixed points of the explicit Euler method. 
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scheme 

c, 

c ; 

E-0 ’ 

0.65 

0.9 

minmod 

0.5 

0.75 

superbee 

— 

0.7 

van Leer 

0.6 

0.7 

van Albada 

0.6 

0.7 


Table 4.1. Convergence regions for the explicit schemes. 


scheme 

converges to x\ 

converges to x-i 

E-0 

c < 11 

c > 22.5 

minmod 

c < 11 

c> 21.5 

superbee 

c < 11 

— 

van Leer 

— 

— 

van Albada 

c < 11 

c > 22.5 


Table 4.2. Convergence regions for the implicit schemes. 



at stable shock 

at unstable shock 

scheme 

X 

Y 

stable 

X 


stable 

explicit E-0 

0.36 

-0.47 

A t < 0.057 


BEH 

— 





mm 

At < 0.1043 





irel 

-0.53 

At < 0.1031 

implicit E-0 

0.37 

-0.48 

VA< 

0.40 

-0.32 

At > 1.0585 




0.52 

0.26 

WAt 





-0.29 

-0.53 

VAt 


Table 4.3. Analytical fixed points of the reduced system. 





















Case 

Grid 

N 

Temporal Behavior 


L 

5 

diverge 

1506 

L 

6 

unsteady chaotic 


L 

7 

unsteady chaotic 


L 

8 

steady oscillatory 


L 

9 

steady mono tonic 

M505 

M 

5 

diverge 

M506 

M 

6 

steady monotonic 

M507 

M 

7 

steady raonotonic 

2T505 

H 

5 

steady monotonic 

H 507 

H 

7 

steady raonotonic 


Case 

Grid 

N 

Temporal Behavior 

B i i M 

mm 

5 

diverge 

■ 

mm 

6 

time periodic 

1 

mm 

7 

time periodic 


wm 

8 

steady oscillatory 


19 

9 

steady raonotonic 

M505 

M 

5 

diverge 

AT506 

M 

6 

steady monotonic 

M507 

M 

7 

steady monotonic 

wjjky| 

H 

5 

steady monotonic 


H 

7 

steady monotonic 


Table 5. 1 . Mesh refinement results for Re — 500 at t < 800. Table 5.2. Mesh refinement results for Re = 500 at t < 2000. 


Case 

Grid 

N 

Temporal Behavior 


L 

9 

diverge 

£811 

L 

11 

unsteady chaotic 


M 

5 

diverge 

M806 

M 

6 

diverge 

M807 

M 

7 

unsteady chaotic 

M808 

M 

8 

steady oscillatory 

M809 

M 

9 

steady monotonic 

if 805 

E 

5 

diverge 

ES 06 

E 

6 

steady raonotonic 

JJ807 

E 

7 

steady raonotonic 

JI808 

E 

8 

steady raonotonic 

J5T809 

E 

9 

steady monotonic 


Table 5.3. Mesh refinement results for Re = 800 at t < 800. 


Case 



Temporal Behavior 


£ 

9 

diverge 

£811 

£ 

11 

time periodic 

Af805 

M 

5 

diverge 

Af806 

M 

6 

diverge 

M807 

M 

7 

diverge 

M808 

M 

8 

steady oscillatory 

AT809 

M 

9 

steady monotonic 

ESOb 

E 

5 

diverge 

J5T806 

E 

6 

steady monotonic 


E 

7 

steady monotonic 


E 

8 

steady raonotonic 


E 

9 

steady raonotonic 


Table 5.4. Mesh refinement results for Re = 800 at t < 2000. 


gg 

Af807 

Af809 

M811 

El 

diverge 

diverge 

diverge 

E 

diverge 

steady raonotonic 

steady monotonic 


diverge 

diverge 

diverge 

m 

diverge 

steady monotonic 

steady raonotonic 

a 

diverge 

diverge 

diverge 


Table 5.5. Initial data study for Re = 800 and M grid at t < 2000. 
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Fig. 3.5a. Stability of steady-state solutions arising through three types of bifurcation phenomena 
( — stable, — unstable). 
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Fig. 3.8. “Full” Bi: 
au(l — u)(6 — u). 




Bifurcation Diagrams & Basins of Attraction 

u’ = au (1-u) 


Modified Euler 


Improved Euler 







Kutta 


R-K 4 



Fig. 3.9. 
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Bifurcation Diagrams & Basins of Attraction 
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Fig. 3.10. 







Bifurcation Diagrams & Basins of Attraction 

u’ = a u (1-u) 
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Fig. 3.12. 




Bifurcation Diagrams & Basins of Attraction 

u' = a u (1-u) 
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Fig. 3.13. 




Bifurcation Diagrams & Basins of Attraction 

u’ = au(1-u) 
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Fig. 3. 14 











Bifurcation Diagrams & Basins of Attraction 
u’ = a u (1 - u) (0.2 - u) 
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Fig. 3.18. 







L2 norm of U Spectral radius 



P2 = OtAt 


Fig. 3.20. Comparison for discretization UP1UI/EE and pi = 7 of (a) linearized stability analysis 
(around u* = 2) and (b) nonlinear solution behavior with initial data U° 

(• : steady-state solution and other asymptotic solution, blank space : divergent 





L2 norm of U Spectral radius 



Fig. 3.21. Comparison for discretization UP1UI/ME and p, = 7 of (a) linearized stability anal- 
ysis (around u * = 2) and (b) nonlinear solution behavior with initial data U° (• 
steady-state solution and other asymptotic solution, blank space : divergent solution) 


1 






PI 

Fig. 3.22. Average wave speed W versus pi of the numerical solution with explicit Euler time 
discretization and spatial discretization UP1PW. 
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CHA/EE 


CHA/ME 



Fig. 3.23. Numerical basin of attraction for discretizations (a) CHA/EE, (b) CHA/ME, (c) 
CHA/LIE and (d) CHA/LT with J = 4,pi = 0.1 and pa = 0.5 (A : exact steady-state 
solution, ■ : spurious steady-state solution, • : other spurious asymptotic solution, 
blank space: divergent solution). 
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CHA/LIE CHA/LT 



Fig. 3.24. Numerical basin of attraction for discretizations (a) CHA/LIE and (b) CHA/LT with 
J = 4, pi = 0.1 and pa = 3 (CFL = 0.3) (A : exact steady-state solution, ■ : spurious 
steady-state solution, • : other spurious asymptotic solution, blank space: divergent 
solution). 


CHA/LIE CHA/LT 



Fig 3 25. Numerical basin of attraction for discretizations (a) CHA/LIE and (b) CHA/LT with 
7 = 4, px =0.1 and pa = 10 (CFL = 1) (A : exact steady-state solution, ■ : spurious 
steady-state solution, • : other spurious asymptotic solution, blank space: divergent 
solution). 
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Log 10< Error In 8) 


Error in entropy Convergence exponent 



Fig. 4.5. Error in entropy (left) and convergence exponent (right) of (4.7) for a second-order 
UNO scheme. 


Ull 


• variable 
■ fixed 



Fig. 4.6. Grid points of the reduced Embid et al. problem. 
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Mean Velocity Profile 


O 



Fig. 5.4. Near wall mean-velocity profiles, o; "corrected" data of Eckelmann (1974); 

lower wall; : upper wall; ; "law of the wall ; u = y , u — 2.5 ln(y +5.5 



Fig. 5.5. Root-mean-square velocity fluctuations normalized by wall shear velocity. 

„ . v • w Cross-channel coordinate normalized by 6 = 2L. 

u nns’ ’ mis’ * nns- 


135 






Fig. 5.6. Reynolds shear stress normalized by wall shear velocity. , -«v; 

-uv ♦ (1 /Re)du/dy, , Total shear stress for fully developed channel 


136 






137 


Fig. 5.9. Streamlines for 1509 (steady) and 1507 (spurious time-periodic) cases (vertical 
scale expanded four times for viewing). 
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Fig. 5. 1 3. Vertical velocity time histories for 1508 and 1509 with At = 0.05 and 0.10. 
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Fig. 5.16. Circular cylinder vertical motion time histories for the lagged structures for A t = 
0.01, 0.02, 0.03, 0.04, 0.05, 0.06 (Re = 500, = 0.2). 
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At-0.01, Lagged Coupling At* 0.0 5, Lagged Co upHng 



Fig. 5.18. Circular cylinder vertical motion time histories for fully coupled and lagged structures for 
At = 0.01 and 0.05. 
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Fig. 5.19. Delta wing configuration. 




Fig. 5.20. Contours of total pressure on vertical plane through vortex center. 
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Fie. 5.21. Lift coefficient time histories computed Fig. 5.22. Frequency spectra of velocity fluctuations 


on Grid 1 for a = 30° using three differ- 
ent initial angles of attack. 


at point within vortex breakdown region 
for ag = 25°. 



for a 0 = 17°. 


Grid 1 



Fig. 5.24. Axial vorticity profile through the vortex 
core upstream of breakdown for the three 
initial angles of attack. 
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Fig. 5.25. Lift coefficient time histories computed 
on Grid 2 for a = 30° using three differ- 
ent initial angles of attack. 
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